ANCOVA in Writing (TDE) (Writing (TDE))

Geiser C. Challco

NOTE:

Setting Initial Variables

dv = "score.tde"
dv.pos = "score.tde.pos"
dv.pre = "score.tde.pre"

fatores2 <- c("Sexo","Zona","Cor.Raca","Serie","score.tde.quintile")
lfatores2 <- as.list(fatores2)
names(lfatores2) <- fatores2

fatores1 <- c("grupo", fatores2)
lfatores1 <- as.list(fatores1)
names(lfatores1) <- fatores1

lfatores <- c(lfatores1)

color <- list()
color[["prepost"]] = c("#ffee65","#f28e2B")
color[["grupo"]] = c("#bcbd22","#008000")
color[["Sexo"]] = c("#FF007F","#4D4DFF")
color[["Zona"]] = c("#AA00FF","#00CCCC")
color[["Cor.Raca"]] = c(
  "Parda"="#b97100","Indígena"="#9F262F",
  "Branca"="#87c498", "Preta"="#848283","Amarela"="#D6B91C"
)

level <- list()
level[["grupo"]] = c("Controle","Experimental")
level[["Sexo"]] = c("F","M")
level[["Zona"]] = c("Rural","Urbana")
level[["Cor.Raca"]] = c("Parda","Indígena","Branca", "Preta","Amarela")
level[["Serie"]] = c("6 ano","7 ano","8 ano","9 ano")

# ..

ymin <- 0
ymax <- 0

ymin.ci <- 0
ymax.ci <- 0


color[["grupo:Sexo"]] = c(
  "Controle:F"="#ff99cb", "Controle:M"="#b7b7ff",
  "Experimental:F"="#FF007F", "Experimental:M"="#4D4DFF",
  "Controle.F"="#ff99cb", "Controle.M"="#b7b7ff",
  "Experimental.F"="#FF007F", "Experimental.M"="#4D4DFF"
)
color[["grupo:Zona"]] = c(
  "Controle:Rural"="#b2efef","Controle:Urbana"="#e5b2ff",
  "Experimental:Rural"="#00CCCC", "Experimental:Urbana"="#AA00FF",
  "Controle.Rural"="#b2efef","Controle.Urbana"="#e5b2ff",
  "Experimental.Rural"="#00CCCC", "Experimental.Urbana"="#AA00FF"
)
color[["grupo:Cor.Raca"]] = c(
    "Controle:Parda"="#e3c699", "Experimental:Parda"="#b97100",
    "Controle:Indígena"="#e2bdc0", "Experimental:Indígena"="#9F262F",
    "Controle:Branca"="#c0e8cb", "Experimental:Branca"="#87c498",
    "Controle:Preta"="#dad9d9", "Experimental:Preta"="#848283",
    "Controle:Amarela"="#eee3a4", "Experimental:Amarela"="#D6B91C",
    
    "Controle.Parda"="#e3c699", "Experimental.Parda"="#b97100",
    "Controle.Indígena"="#e2bdc0", "Experimental.Indígena"="#9F262F",
    "Controle.Branca"="#c0e8cb", "Experimental.Branca"="#87c498",
    "Controle.Preta"="#dad9d9", "Experimental.Preta"="#848283",
    "Controle.Amarela"="#eee3a4", "Experimental.Amarela"="#D6B91C"
)


for (coln in c("vocab","vocab.teach","vocab.non.teach","score.tde",
               "TFL.lidas.per.min","TFL.corretas.per.min","TFL.erradas.per.min","TFL.omitidas.per.min",
               "leitura.compreensao")) {
  color[[paste0(coln,".quintile")]] = c("#BF0040","#FF0000","#800080","#0000FF","#4000BF")
  level[[paste0(coln,".quintile")]] = c("1st quintile","2nd quintile","3rd quintile","4th quintile","5th quintile")
  color[[paste0("grupo:",coln,".quintile")]] = c(
    "Experimental.1st quintile"="#BF0040", "Controle.1st quintile"="#d8668c",
    "Experimental.2nd quintile"="#FF0000", "Controle.2nd quintile"="#ff7f7f",
    "Experimental.3rd quintile"="#8fce00", "Controle.3rd quintile"="#ddf0b2",
    "Experimental.4th quintile"="#0000FF", "Controle.4th quintile"="#b2b2ff",
    "Experimental.5th quintile"="#4000BF", "Controle.5th quintile"="#b299e5",
    
    "Experimental:1st quintile"="#BF0040", "Controle:1st quintile"="#d8668c",
    "Experimental:2nd quintile"="#FF0000", "Controle:2nd quintile"="#ff7f7f",
    "Experimental:3rd quintile"="#8fce00", "Controle:3rd quintile"="#ddf0b2",
    "Experimental:4th quintile"="#0000FF", "Controle:4th quintile"="#b2b2ff",
    "Experimental:5th quintile"="#4000BF", "Controle:5th quintile"="#b299e5")
}


tdat <- read_excel("../data/data.xlsx", sheet = "sumary")
tdat <- tdat[!is.na(tdat[["WG.Grupo"]]),]
tdat$grupo <- factor(tdat[["WG.Grupo"]], level[["grupo"]]) 

gdat <- tdat[which(is.na(tdat$Necessidade.Deficiencia) & !is.na(tdat$WG.Grupo)),]



dat <- gdat
dat$grupo <- factor(dat[["WG.Grupo"]], level[["grupo"]])
for (coln in c(names(lfatores))) {
  dat[[coln]] <- factor(dat[[coln]], level[[coln]][level[[coln]] %in% unique(dat[[coln]])])
}
dat <- dat[which(!is.na(dat[[dv.pre]]) & !is.na(dat[[dv.pos]])),]
dat <- dat[,c("id",names(lfatores),dv.pre,dv.pos)]

dat.long <- rbind(dat, dat)
dat.long$time <- c(rep("pre", nrow(dat)), rep("pos", nrow(dat)))
dat.long$time <- factor(dat.long$time, c("pre","pos"))
dat.long[[dv]] <- c(dat[[dv.pre]], dat[[dv.pos]])


for (f in c("grupo", names(lfatores))) {
  if (is.null(color[[f]]) && length(unique(dat[[f]])) > 0) 
      color[[f]] <- distinctColorPalette(length(unique(dat[[f]])))
}
for (f in c(fatores2)) {
  if (is.null(color[[paste0("grupo:",f)]]) && length(unique(dat[[f]])) > 0)
    color[[paste0("grupo:",f)]] <- distinctColorPalette(length(unique(dat[["grupo"]]))*length(unique(dat[[f]])))
}

ldat <- list()
laov <- list()
lpwc <- list()
lemms <- list()

Descriptive Statistics of Initial Data

df <- get.descriptives(dat, c(dv.pre, dv.pos), c("grupo"), 
                       include.global = T, symmetry.test = T, normality.test = F)
df <- plyr::rbind.fill(
  df, do.call(plyr::rbind.fill, lapply(lfatores2, FUN = function(f) {
    if (nrow(dat) > 0 && sum(!is.na(unique(dat[[f]]))) > 1)
      get.descriptives(dat, c(dv.pre,dv.pos), c("grupo", f),
                       symmetry.test = T, normality.test = F)
    }))
)
## Warning: There was 1 warning in `mutate()`.
## ℹ In argument: `ci = abs(stats::qt(alpha/2, .data$n - 1) * .data$se)`.
## Caused by warning:
## ! There was 1 warning in `mutate()`.
## ℹ In argument: `ci = abs(stats::qt(alpha/2, .data$n - 1) * .data$se)`.
## Caused by warning in `stats::qt()`:
## ! NaNs produced
## There was 1 warning in `mutate()`.
## ℹ In argument: `ci = abs(stats::qt(alpha/2, .data$n - 1) * .data$se)`.
## Caused by warning:
## ! There was 1 warning in `mutate()`.
## ℹ In argument: `ci = abs(stats::qt(alpha/2, .data$n - 1) * .data$se)`.
## Caused by warning in `stats::qt()`:
## ! NaNs produced
df <- df[,c(fatores1[fatores1 %in% colnames(df)],"variable",
            colnames(df)[!colnames(df) %in% c(fatores1,"variable")])]
grupo Sexo Zona Cor.Raca Serie score.tde.quintile variable n mean median min max sd se ci iqr symmetry skewness kurtosis
Controle score.tde.pre 485 37.631 44.0 0 75 19.176 0.871 1.711 33.00 YES -0.443 -1.017
Experimental score.tde.pre 636 36.748 39.0 0 73 16.891 0.670 1.315 24.00 YES -0.303 -0.711
score.tde.pre 1121 37.130 41.0 0 75 17.912 0.535 1.050 28.00 YES -0.370 -0.857
Controle score.tde.pos 485 33.816 36.0 0 73 20.988 0.953 1.873 37.00 YES -0.108 -1.287
Experimental score.tde.pos 636 35.376 38.0 0 74 18.900 0.749 1.472 30.00 YES -0.158 -1.024
score.tde.pos 1121 34.701 37.0 0 74 19.836 0.592 1.162 33.00 YES -0.145 -1.138
Controle F score.tde.pre 247 40.652 46.0 0 72 17.259 1.098 2.163 24.50 NO -0.692 -0.519
Controle M score.tde.pre 238 34.496 42.0 0 75 20.552 1.332 2.624 37.00 YES -0.173 -1.308
Experimental F score.tde.pre 319 39.408 43.0 0 73 16.802 0.941 1.851 25.00 YES -0.460 -0.599
Experimental M score.tde.pre 317 34.073 36.0 0 71 16.580 0.931 1.832 23.00 YES -0.175 -0.723
Controle F score.tde.pos 247 36.567 42.0 0 72 19.883 1.265 2.492 35.00 YES -0.290 -1.208
Controle M score.tde.pos 238 30.962 28.5 0 73 21.752 1.410 2.778 37.00 YES 0.098 -1.291
Experimental F score.tde.pos 319 38.138 42.0 0 73 19.397 1.086 2.137 32.00 YES -0.311 -1.059
Experimental M score.tde.pos 317 32.596 34.0 0 74 17.992 1.011 1.988 29.00 YES -0.050 -0.909
Controle Rural score.tde.pre 243 36.181 43.0 0 69 19.626 1.259 2.480 34.50 YES -0.429 -1.220
Controle Urbana score.tde.pre 109 39.468 43.0 1 75 18.034 1.727 3.424 26.00 YES -0.361 -0.676
Controle score.tde.pre 133 38.774 45.0 0 73 19.183 1.663 3.290 32.00 YES -0.480 -1.003
Experimental Rural score.tde.pre 284 34.592 35.5 0 73 17.467 1.036 2.040 25.25 YES -0.140 -0.823
Experimental Urbana score.tde.pre 167 37.850 42.0 0 71 16.633 1.287 2.541 22.00 YES -0.491 -0.519
Experimental score.tde.pre 185 39.065 42.0 0 69 15.877 1.167 2.303 23.00 YES -0.345 -0.691
Controle Rural score.tde.pos 243 34.243 39.0 0 72 21.515 1.380 2.719 38.00 YES -0.201 -1.321
Controle Urbana score.tde.pos 109 29.835 26.0 0 73 20.374 1.951 3.868 37.00 YES 0.237 -1.198
Controle score.tde.pos 133 36.301 39.0 0 73 20.179 1.750 3.461 34.00 YES -0.200 -1.185
Experimental Rural score.tde.pos 284 33.627 33.0 0 72 19.232 1.141 2.246 33.00 YES 0.039 -1.088
Experimental Urbana score.tde.pos 167 34.108 38.0 0 74 18.899 1.462 2.887 30.50 YES -0.228 -1.049
Experimental score.tde.pos 185 39.205 42.0 0 73 17.909 1.317 2.598 26.00 YES -0.387 -0.775
Controle Parda score.tde.pre 162 36.741 43.0 0 66 18.900 1.485 2.932 31.25 NO -0.595 -0.975
Controle Indígena score.tde.pre 11 42.000 46.0 4 65 17.297 5.215 11.621 9.50 NO -1.005 -0.104
Controle Branca score.tde.pre 50 41.460 46.0 2 67 16.942 2.396 4.815 14.50 NO -0.821 -0.311
Controle score.tde.pre 262 37.267 43.0 0 75 19.798 1.223 2.408 34.00 YES -0.268 -1.181
Experimental Parda score.tde.pre 190 34.253 35.0 0 68 17.537 1.272 2.510 27.00 YES -0.133 -0.931
Experimental Indígena score.tde.pre 15 27.067 26.0 6 54 17.065 4.406 9.450 27.00 YES 0.126 -1.422
Experimental Branca score.tde.pre 61 34.623 36.0 0 73 19.728 2.526 5.053 33.00 YES -0.045 -1.037
Experimental Amarela score.tde.pre 1 23.000 23.0 23 23 0.00 few data 0.000 0.000
Experimental score.tde.pre 369 38.816 42.0 0 71 15.723 0.819 1.610 21.00 YES -0.438 -0.449
Controle Parda score.tde.pos 162 32.463 35.0 0 68 20.798 1.634 3.227 39.75 YES -0.166 -1.437
Controle Indígena score.tde.pos 11 45.636 52.0 4 63 19.931 6.009 13.390 19.50 NO -0.968 -0.645
Controle Branca score.tde.pos 50 34.300 38.0 0 72 20.518 2.902 5.831 34.75 YES -0.208 -1.184
Controle score.tde.pos 262 34.065 34.0 0 73 21.180 1.309 2.577 36.00 YES -0.024 -1.252
Experimental Parda score.tde.pos 190 32.447 30.0 0 71 19.467 1.412 2.786 32.75 YES 0.082 -1.112
Experimental Indígena score.tde.pos 15 30.333 26.0 1 57 17.715 4.574 9.810 30.50 YES 0.117 -1.542
Experimental Branca score.tde.pos 61 33.311 34.0 0 72 21.189 2.713 5.427 34.00 YES 0.076 -1.201
Experimental Amarela score.tde.pos 1 23.000 23.0 23 23 0.00 few data 0.000 0.000
Experimental score.tde.pos 369 37.463 40.0 0 74 18.045 0.939 1.847 27.00 YES -0.336 -0.844
Controle 6 ano score.tde.pre 134 28.993 26.5 0 63 19.526 1.687 3.336 36.00 YES 0.057 -1.526
Controle 7 ano score.tde.pre 141 34.716 41.0 0 69 17.842 1.503 2.971 29.00 YES -0.265 -1.162
Controle 8 ano score.tde.pre 89 42.584 46.0 0 72 16.382 1.737 3.451 17.00 NO -0.833 0.214
Controle 9 ano score.tde.pre 121 46.950 51.0 0 75 17.124 1.557 3.082 15.00 NO -1.124 0.589
Experimental 6 ano score.tde.pre 159 29.000 29.0 0 65 15.950 1.265 2.498 24.50 YES -0.043 -0.909
Experimental 7 ano score.tde.pre 187 36.540 38.0 0 68 14.891 1.089 2.148 21.50 YES -0.195 -0.682
Experimental 8 ano score.tde.pre 143 39.000 44.0 0 73 18.129 1.516 2.997 25.50 NO -0.515 -0.734
Experimental 9 ano score.tde.pre 147 43.204 45.0 0 71 15.845 1.307 2.583 22.00 NO -0.641 -0.044
Controle 6 ano score.tde.pos 134 24.231 19.5 0 66 20.155 1.741 3.444 40.00 YES 0.381 -1.295
Controle 7 ano score.tde.pos 141 28.362 22.0 0 72 20.646 1.739 3.437 37.00 YES 0.379 -1.190
Controle 8 ano score.tde.pos 89 39.787 43.0 0 71 17.900 1.897 3.771 28.00 NO -0.521 -0.552
Controle 9 ano score.tde.pos 121 46.397 50.0 0 73 16.448 1.495 2.960 20.00 NO -0.786 0.049
Experimental 6 ano score.tde.pos 159 24.094 22.0 0 65 16.759 1.329 2.625 26.00 NO 0.518 -0.626
Experimental 7 ano score.tde.pos 187 33.524 36.0 0 69 17.147 1.254 2.474 27.50 YES -0.063 -0.934
Experimental 8 ano score.tde.pos 143 41.042 45.0 0 73 19.177 1.604 3.170 26.50 NO -0.568 -0.652
Experimental 9 ano score.tde.pos 147 44.422 48.0 0 74 16.143 1.331 2.631 22.50 NO -0.733 0.081
Controle 1st quintile score.tde.pre 113 9.230 10.0 0 18 5.560 0.523 1.036 9.00 YES -0.093 -1.227
Controle 2nd quintile score.tde.pre 59 24.542 25.0 19 31 3.780 0.492 0.985 6.00 YES -0.056 -1.175
Controle 3rd quintile score.tde.pre 54 38.815 39.5 32 42 3.004 0.409 0.820 5.00 NO -0.510 -1.049
Controle 4th quintile score.tde.pre 135 47.370 47.0 43 51 2.579 0.222 0.439 4.50 YES -0.238 -1.107
Controle 5th quintile score.tde.pre 124 58.621 58.0 52 75 5.367 0.482 0.954 8.00 NO 0.794 0.036
Experimental 1st quintile score.tde.pre 112 9.777 10.0 0 18 5.445 0.515 1.020 9.25 YES -0.202 -1.150
Experimental 2nd quintile score.tde.pre 117 25.855 26.0 19 31 3.422 0.316 0.627 6.00 YES -0.218 -1.044
Experimental 3rd quintile score.tde.pre 141 37.383 38.0 32 42 3.177 0.268 0.529 6.00 YES -0.147 -1.221
Experimental 4th quintile score.tde.pre 132 46.697 47.0 43 51 2.489 0.217 0.429 4.25 YES 0.144 -1.180
Experimental 5th quintile score.tde.pre 134 58.336 57.5 52 73 5.029 0.434 0.859 8.00 NO 0.688 -0.352
Controle 1st quintile score.tde.pos 113 8.522 5.0 0 44 8.588 0.808 1.601 12.00 NO 1.083 1.171
Controle 2nd quintile score.tde.pos 59 19.966 20.0 0 66 12.947 1.686 3.374 13.50 NO 1.282 2.560
Controle 3rd quintile score.tde.pos 54 30.185 28.5 2 62 15.298 2.082 4.176 22.25 YES 0.232 -0.869
Controle 4th quintile score.tde.pos 135 43.215 46.0 14 63 11.360 0.978 1.934 16.00 NO -0.680 -0.300
Controle 5th quintile score.tde.pos 124 54.806 56.0 22 73 10.905 0.979 1.938 11.25 NO -0.916 0.870
Experimental 1st quintile score.tde.pos 112 11.384 10.0 0 51 9.500 0.898 1.779 12.00 NO 1.241 2.205
Experimental 2nd quintile score.tde.pos 117 24.675 25.0 0 60 12.203 1.128 2.234 16.00 YES 0.028 -0.429
Experimental 3rd quintile score.tde.pos 141 34.872 37.0 0 64 13.013 1.096 2.167 18.00 YES -0.298 -0.483
Experimental 4th quintile score.tde.pos 132 45.030 47.0 10 69 10.621 0.924 1.829 11.00 NO -0.801 0.896
Experimental 5th quintile score.tde.pos 134 55.791 57.0 14 74 10.619 0.917 1.814 12.00 NO -1.174 2.100

ANCOVA and Pairwise for one factor: grupo

Without remove non-normal data

pdat = remove_group_data(dat[!is.na(dat[["grupo"]]),], "score.tde.pos", "grupo")

pdat.long <- rbind(pdat[,c("id","grupo")], pdat[,c("id","grupo")])
pdat.long[["time"]] <- c(rep("pre", nrow(pdat)), rep("pos", nrow(pdat)))
pdat.long[["time"]] <- factor(pdat.long[["time"]], c("pre","pos"))
pdat.long[["score.tde"]] <- c(pdat[["score.tde.pre"]], pdat[["score.tde.pos"]])

aov = anova_test(pdat, score.tde.pos ~ score.tde.pre + grupo)
laov[["grupo"]] <- get_anova_table(aov)
pwc <- emmeans_test(pdat, score.tde.pos ~ grupo, covariate = score.tde.pre,
                    p.adjust.method = "bonferroni")
pwc.long <- emmeans_test(dplyr::group_by_at(pdat.long, "grupo"),
                          score.tde ~ time,
                          p.adjust.method = "bonferroni")
lpwc[["grupo"]] <- plyr::rbind.fill(pwc, pwc.long)
ds <- get.descriptives(pdat, "score.tde.pos", "grupo", covar = "score.tde.pre")
ds <- merge(ds[ds$variable != "score.tde.pre",],
            ds[ds$variable == "score.tde.pre", !colnames(ds) %in% c("variable")],
            by = "grupo", all.x = T, suffixes = c("", ".score.tde.pre"))
ds <- merge(get_emmeans(pwc), ds, by = "grupo", suffixes = c(".emms", ""))
ds <- ds[,c("grupo","n","mean.score.tde.pre","se.score.tde.pre","mean","se",
            "emmean","se.emms","conf.low","conf.high")]

colnames(ds) <- c("grupo", "N", paste0(c("M","SE")," (pre)"),
                  paste0(c("M","SE"), " (unadj)"),
                  paste0(c("M", "SE"), " (adj)"), "conf.low", "conf.high")

lemms[["grupo"]] <- ds

Computing ANCOVA and PairWise After removing non-normal data (OK)

wdat = pdat 

res = residuals(lm(score.tde.pos ~ score.tde.pre + grupo, data = wdat))
non.normal = getNonNormal(res, wdat$id, plimit = 0.05)

wdat = wdat[!wdat$id %in% non.normal,]

wdat.long <- rbind(wdat[,c("id","grupo")], wdat[,c("id","grupo")])
wdat.long[["time"]] <- c(rep("pre", nrow(wdat)), rep("pos", nrow(wdat)))
wdat.long[["time"]] <- factor(wdat.long[["time"]], c("pre","pos"))
wdat.long[["score.tde"]] <- c(wdat[["score.tde.pre"]], wdat[["score.tde.pos"]])

ldat[["grupo"]] = wdat

(non.normal)
##   [1] "P962"  "P1128" "P984"  "P1129" "P3637" "P1126" "P1971" "P572"  "P3021"
##  [10] "P2871" "P1139" "P921"  "P2858" "P3054" "P3019" "P2995" "P2848" "P2870"
##  [19] "P1117" "P2929" "P929"  "P2835" "P2983" "P2946" "P2865" "P2910" "P3007"
##  [28] "P3015" "P2964" "P2975" "P2876" "P2937" "P2831" "P2953" "P2986" "P1018"
##  [37] "P2950" "P2883" "P1878" "P2917" "P2866" "P2997" "P2880" "P908"  "P2978"
##  [46] "P2864" "P3005" "P2974" "P3026" "P2904" "P2994" "P1111" "P2846" "P2913"
##  [55] "P2886" "P2905" "P3533" "P914"  "P2888" "P2867" "P1983" "P3574" "P1840"
##  [64] "P3545" "P2993" "P976"  "P2973" "P2967" "P3029" "P3548" "P3000" "P2959"
##  [73] "P2956" "P2982" "P3475" "P2854" "P1068" "P3476" "P3674" "P3660" "P2004"
##  [82] "P3666" "P2190" "P2839" "P2861" "P1885" "P2947" "P606"  "P1149" "P609" 
##  [91] "P2843" "P3228" "P2845" "P2868" "P2860" "P3024" "P2852" "P2990" "P2969"
## [100] "P809"  "P2909" "P924"  "P1914" "P1056" "P2847" "P1046" "P1118" "P541" 
## [109] "P3651" "P3221" "P2903" "P2979" "P3537" "P2981" "P3237" "P1127" "P1916"
## [118] "P2832" "P1828" "P3571" "P996"  "P2869" "P2836" "P2879" "P1900" "P2891"
## [127] "P3608" "P2882" "P2951" "P1673" "P993"  "P1910" "P548"  "P492"  "P1887"
## [136] "P848"  "P1964" "P859"  "P2971" "P2922"
aov = anova_test(wdat, score.tde.pos ~ score.tde.pre + grupo)
laov[["grupo"]] <- merge(get_anova_table(aov), laov[["grupo"]],
                            by="Effect", suffixes = c("","'"))

(df = get_anova_table(aov))
## ANOVA Table (type II tests)
## 
##          Effect DFn DFd        F        p p<.05   ges
## 1 score.tde.pre   1 978 5874.442 0.00e+00     * 0.857
## 2         grupo   1 978   24.670 8.03e-07     * 0.025
Effect DFn DFd F p p<.05 ges
score.tde.pre 1 978 5874.442 0 * 0.857
grupo 1 978 24.670 0 * 0.025 
pwc <- emmeans_test(wdat, score.tde.pos ~ grupo, covariate = score.tde.pre,
                    p.adjust.method = "bonferroni")
term .y. group1 group2 df statistic p p.adj p.adj.signif
score.tde.pre*grupo score.tde.pos Controle Experimental 978 -4.967 0 0 **** 
pwc.long <- emmeans_test(dplyr::group_by_at(wdat.long, "grupo"),
                         score.tde ~ time,
                         p.adjust.method = "bonferroni")
lpwc[["grupo"]] <- merge(plyr::rbind.fill(pwc, pwc.long), lpwc[["grupo"]],
                            by=c("grupo","term",".y.","group1","group2"),
                            suffixes = c("","'"))
grupo term .y. group1 group2 df statistic p p.adj p.adj.signif
Controle time score.tde pre pos 1958 1.465 0.143 0.143 ns
Experimental time score.tde pre pos 1958 -0.402 0.688 0.688 ns 
ds <- get.descriptives(wdat, "score.tde.pos", "grupo", covar = "score.tde.pre")
ds <- merge(ds[ds$variable != "score.tde.pre",],
            ds[ds$variable == "score.tde.pre", !colnames(ds) %in% c("variable")],
            by = "grupo", all.x = T, suffixes = c("", ".score.tde.pre"))
ds <- merge(get_emmeans(pwc), ds, by = "grupo", suffixes = c(".emms", ""))
ds <- ds[,c("grupo","n","mean.score.tde.pre","se.score.tde.pre","mean","se",
            "emmean","se.emms","conf.low","conf.high")]

colnames(ds) <- c("grupo", "N", paste0(c("M","SE")," (pre)"),
                  paste0(c("M","SE"), " (unadj)"),
                  paste0(c("M", "SE"), " (adj)"), "conf.low", "conf.high")

lemms[["grupo"]] <- merge(ds, lemms[["grupo"]], by=c("grupo"), suffixes = c("","'"))
grupo N M (pre) SE (pre) M (unadj) SE (unadj) M (adj) SE (adj) conf.low conf.high
Controle 424 37.351 0.975 35.432 1.018 35.112 0.358 34.409 35.815
Experimental 557 36.770 0.741 37.230 0.777 37.473 0.313 36.860 38.087

Plots for ancova

plots <- oneWayAncovaPlots(
  wdat, "score.tde.pos", "grupo", aov, list("grupo"=pwc), addParam = c("mean_ci"),
  font.label.size=10, step.increase=0.05, p.label="p.adj",
  subtitle = which(aov$Effect == "grupo"))
if (!is.null(nrow(plots[["grupo"]]$data)))
  plots[["grupo"]]   +
  ggplot2::ylab("Writing (TDE)") +
  theme(axis.title = element_text(size = 14),
        legend.text = element_text(size = 16),
        plot.subtitle = element_text(size = 18)) +
  if (ymin.ci < ymax.ci) ggplot2::ylim(ymin.ci, ymax.ci)

plots <- oneWayAncovaBoxPlots(
  wdat, "score.tde.pos", "grupo", aov, pwc, covar = "score.tde.pre",
  theme = "classic", color = color[["grupo"]],
  subtitle = which(aov$Effect == "grupo"))
if (length(unique(wdat[["grupo"]])) > 1)
  plots[["grupo"]] + ggplot2::ylab("Writing (TDE)") +
  ggplot2::scale_x_discrete(labels=c('pre', 'pos')) +
  if (ymin < ymax) ggplot2::ylim(ymin, ymax)

if (length(unique(wdat.long[["grupo"]])) > 1)
  plots <- oneWayAncovaBoxPlots(
    wdat.long, "score.tde", "grupo", aov, pwc.long,
    pre.post = "time", theme = "classic", color = color$prepost)
if (length(unique(wdat.long[["grupo"]])) > 1)
  plots[["grupo"]] + ggplot2::ylab("Writing (TDE)") +
  if (ymin < ymax) ggplot2::ylim(ymin, ymax)

Checking linearity assumption

ggscatter(wdat, x = "score.tde.pre", y = "score.tde.pos", size = 0.5,
          color = "grupo", add = "reg.line")+
  stat_regline_equation(
    aes(label =  paste(..eq.label.., ..rr.label.., sep = "~~~~"), color = grupo)
  ) +
  ggplot2::labs(subtitle = rstatix::get_test_label(aov, detailed = T, row = which(aov$Effect == "grupo"))) +
  ggplot2::scale_color_manual(values = color[["grupo"]]) +
  ggplot2::xlab("Writing (TDE)") +
  ggplot2::ylab("Writing (TDE)") +
  theme(axis.title = element_text(size = 14),
        legend.text = element_text(size = 16),
        plot.subtitle = element_text(size = 18)) +
  if (ymin < ymax) ggplot2::ylim(ymin, ymax)
## Warning: The dot-dot notation (`..eq.label..`) was deprecated in ggplot2 3.4.0.
## ℹ Please use `after_stat(eq.label)` instead.
## ℹ The deprecated feature was likely used in the ggpubr package.
##   Please report the issue at <https://github.com/kassambara/ggpubr/issues>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

Checking normality and homogeneity

res <- augment(lm(score.tde.pos ~ score.tde.pre + grupo, data = wdat))
shapiro_test(res$.resid)
## # A tibble: 1 × 3
##   variable   statistic p.value
##   <chr>          <dbl>   <dbl>
## 1 res$.resid     0.996  0.0185
levene_test(res, .resid ~ grupo)
## # A tibble: 1 × 4
##     df1   df2 statistic     p
##   <int> <int>     <dbl> <dbl>
## 1     1   979   0.00118 0.973

ANCOVA and Pairwise for two factors grupo:Sexo

Without remove non-normal data

pdat = remove_group_data(dat[!is.na(dat[["grupo"]]) & !is.na(dat[["Sexo"]]),],
                         "score.tde.pos", c("grupo","Sexo"))
pdat = pdat[pdat[["Sexo"]] %in% do.call(
  intersect, lapply(unique(pdat[["grupo"]]), FUN = function(x) {
    unique(pdat[["Sexo"]][which(pdat[["grupo"]] == x)])
  })),]
pdat[["grupo"]] = factor(pdat[["grupo"]], level[["grupo"]])
pdat[["Sexo"]] = factor(
  pdat[["Sexo"]],
  level[["Sexo"]][level[["Sexo"]] %in% unique(pdat[["Sexo"]])])

pdat.long <- rbind(pdat[,c("id","grupo","Sexo")], pdat[,c("id","grupo","Sexo")])
pdat.long[["time"]] <- c(rep("pre", nrow(pdat)), rep("pos", nrow(pdat)))
pdat.long[["time"]] <- factor(pdat.long[["time"]], c("pre","pos"))
pdat.long[["score.tde"]] <- c(pdat[["score.tde.pre"]], pdat[["score.tde.pos"]])

if (length(unique(pdat[["Sexo"]])) >= 2) {
  aov = anova_test(pdat, score.tde.pos ~ score.tde.pre + grupo*Sexo)
  laov[["grupo:Sexo"]] <- get_anova_table(aov)
}
if (length(unique(pdat[["Sexo"]])) >= 2) {
  pwcs <- list()
  pwcs[["Sexo"]] <- emmeans_test(
    group_by(pdat, grupo), score.tde.pos ~ Sexo,
    covariate = score.tde.pre, p.adjust.method = "bonferroni")
  pwcs[["grupo"]] <- emmeans_test(
    group_by(pdat, Sexo), score.tde.pos ~ grupo,
    covariate = score.tde.pre, p.adjust.method = "bonferroni")
  
  pwc <- plyr::rbind.fill(pwcs[["grupo"]], pwcs[["Sexo"]])
  pwc <- pwc[,c("grupo","Sexo", colnames(pwc)[!colnames(pwc) %in% c("grupo","Sexo")])]
}
if (length(unique(pdat[["Sexo"]])) >= 2) {
  pwc.long <- emmeans_test(dplyr::group_by_at(pdat.long, c("grupo","Sexo")),
                           score.tde ~ time,
                           p.adjust.method = "bonferroni")
  lpwc[["grupo:Sexo"]] <- plyr::rbind.fill(pwc, pwc.long)
}
if (length(unique(pdat[["Sexo"]])) >= 2) {
  ds <- get.descriptives(pdat, "score.tde.pos", c("grupo","Sexo"), covar = "score.tde.pre")
  ds <- merge(ds[ds$variable != "score.tde.pre",],
              ds[ds$variable == "score.tde.pre", !colnames(ds) %in% c("variable")],
              by = c("grupo","Sexo"), all.x = T, suffixes = c("", ".score.tde.pre"))
  ds <- merge(get_emmeans(pwcs[["grupo"]]), ds,
              by = c("grupo","Sexo"), suffixes = c(".emms", ""))
  ds <- ds[,c("grupo","Sexo","n","mean.score.tde.pre","se.score.tde.pre","mean","se",
              "emmean","se.emms","conf.low","conf.high")]
  
  colnames(ds) <- c("grupo","Sexo", "N", paste0(c("M","SE")," (pre)"),
                    paste0(c("M","SE"), " (unadj)"),
                    paste0(c("M", "SE"), " (adj)"), "conf.low", "conf.high")
  
  lemms[["grupo:Sexo"]] <- ds
}

Computing ANCOVA and PairWise After removing non-normal data (OK)

if (length(unique(pdat[["Sexo"]])) >= 2) {
  wdat = pdat 
  
  res = residuals(lm(score.tde.pos ~ score.tde.pre + grupo*Sexo, data = wdat))
  non.normal = getNonNormal(res, wdat$id, plimit = 0.05)
  
  wdat = wdat[!wdat$id %in% non.normal,]
  
  wdat.long <- rbind(wdat[,c("id","grupo","Sexo")], wdat[,c("id","grupo","Sexo")])
  wdat.long[["time"]] <- c(rep("pre", nrow(wdat)), rep("pos", nrow(wdat)))
  wdat.long[["time"]] <- factor(wdat.long[["time"]], c("pre","pos"))
  wdat.long[["score.tde"]] <- c(wdat[["score.tde.pre"]], wdat[["score.tde.pos"]])
  
  
  ldat[["grupo:Sexo"]] = wdat
  
  (non.normal)
}
##   [1] "P962"  "P1128" "P984"  "P1129" "P3637" "P1126" "P1971" "P2995" "P2871"
##  [10] "P572"  "P921"  "P3019" "P2994" "P2880" "P2848" "P2983" "P908"  "P3021"
##  [19] "P3054" "P1139" "P2858" "P2835" "P2904" "P1117" "P3015" "P2870" "P2929"
##  [28] "P2910" "P2886" "P914"  "P2913" "P1111" "P2974" "P2964" "P2867" "P2876"
##  [37] "P2937" "P1840" "P2973" "P3007" "P2861" "P2969" "P2953" "P2986" "P2883"
##  [46] "P2831" "P2975" "P1018" "P3476" "P2997" "P3005" "P1878" "P2866" "P2950"
##  [55] "P2993" "P3545" "P2917" "P2978" "P3029" "P2864" "P2888" "P976"  "P2854"
##  [64] "P2959" "P3548" "P2190" "P929"  "P3574" "P2947" "P2845" "P1983" "P3674"
##  [73] "P1149" "P1068" "P3475" "P2004" "P2860" "P2982" "P2839" "P3024" "P1885"
##  [82] "P2868" "P609"  "P2990" "P2946" "P1056" "P2933" "P3020" "P2879" "P873" 
##  [91] "P1887" "P2865" "P3026" "P3660" "P2956" "P3080" "P3666" "P2846" "P3533"
## [100] "P3228" "P2948" "P2850" "P3000" "P809"  "P2905" "P1046" "P1118" "P3221"
## [109] "P1665" "P2903" "P1597" "P2836" "P3537" "P584"  "P854"  "P1127" "P1804"
## [118] "P2869" "P2979" "P3571" "P1652" "P924"  "P1803" "P3088" "P1592" "P2901"
## [127] "P3608" "P848"  "P1964" "P2971" "P1900" "P996"  "P3732" "P492"  "P2922"
## [136] "P859"  "P2980" "P1916" "P2832" "P1858" "P1151" "P2375" "P3156" "P3030"
## [145] "P3061" "P2327" "P2894"
if (length(unique(pdat[["Sexo"]])) >= 2) {
  aov = anova_test(wdat, score.tde.pos ~ score.tde.pre + grupo*Sexo)
  laov[["grupo:Sexo"]] <- merge(get_anova_table(aov), laov[["grupo:Sexo"]],
                                         by="Effect", suffixes = c("","'"))
  df = get_anova_table(aov)
}
Effect DFn DFd F p p<.05 ges
score.tde.pre 1 969 5974.401 0.000 * 0.860
grupo 1 969 23.992 0.000 * 0.024
Sexo 1 969 0.058 0.810 0.000
grupo:Sexo 1 969 0.548 0.459 0.001 
if (length(unique(pdat[["Sexo"]])) >= 2) {
  pwcs <- list()
  pwcs[["Sexo"]] <- emmeans_test(
    group_by(wdat, grupo), score.tde.pos ~ Sexo,
    covariate = score.tde.pre, p.adjust.method = "bonferroni")
  pwcs[["grupo"]] <- emmeans_test(
    group_by(wdat, Sexo), score.tde.pos ~ grupo,
    covariate = score.tde.pre, p.adjust.method = "bonferroni")
  
  pwc <- plyr::rbind.fill(pwcs[["grupo"]], pwcs[["Sexo"]])
  pwc <- pwc[,c("grupo","Sexo", colnames(pwc)[!colnames(pwc) %in% c("grupo","Sexo")])]
}
grupo Sexo term .y. group1 group2 df statistic p p.adj p.adj.signif
F score.tde.pre*grupo score.tde.pos Controle Experimental 969 -3.977 0.000 0.000 ****
M score.tde.pre*grupo score.tde.pos Controle Experimental 969 -2.953 0.003 0.003 **
Controle score.tde.pre*Sexo score.tde.pos F M 969 -0.398 0.691 0.691 ns
Experimental score.tde.pre*Sexo score.tde.pos F M 969 0.663 0.507 0.507 ns 
if (length(unique(pdat[["Sexo"]])) >= 2) {
  pwc.long <- emmeans_test(dplyr::group_by_at(wdat.long, c("grupo","Sexo")),
                           score.tde ~ time,
                           p.adjust.method = "bonferroni")
  lpwc[["grupo:Sexo"]] <- merge(plyr::rbind.fill(pwc, pwc.long),
                                         lpwc[["grupo:Sexo"]],
                                         by=c("grupo","Sexo","term",".y.","group1","group2"),
                                         suffixes = c("","'"))
}
grupo Sexo term .y. group1 group2 df statistic p p.adj p.adj.signif
Controle F time score.tde pre pos 1940 1.089 0.276 0.276 ns
Controle M time score.tde pre pos 1940 0.837 0.403 0.403 ns
Experimental F time score.tde pre pos 1940 -0.403 0.687 0.687 ns
Experimental M time score.tde pre pos 1940 -0.251 0.802 0.802 ns 
if (length(unique(pdat[["Sexo"]])) >= 2) {
  ds <- get.descriptives(wdat, "score.tde.pos", c("grupo","Sexo"), covar = "score.tde.pre")
  ds <- merge(ds[ds$variable != "score.tde.pre",],
              ds[ds$variable == "score.tde.pre", !colnames(ds) %in% c("variable")],
              by = c("grupo","Sexo"), all.x = T, suffixes = c("", ".score.tde.pre"))
  ds <- merge(get_emmeans(pwcs[["grupo"]]), ds,
              by = c("grupo","Sexo"), suffixes = c(".emms", ""))
  ds <- ds[,c("grupo","Sexo","n","mean.score.tde.pre","se.score.tde.pre",
              "mean","se","emmean","se.emms","conf.low","conf.high")]
  
  colnames(ds) <- c("grupo","Sexo", "N", paste0(c("M","SE")," (pre)"),
                    paste0(c("M","SE"), " (unadj)"),
                    paste0(c("M", "SE"), " (adj)"), "conf.low", "conf.high")
  
  lemms[["grupo:Sexo"]] <- merge(ds, lemms[["grupo:Sexo"]],
                                          by=c("grupo","Sexo"), suffixes = c("","'"))
}
grupo Sexo N M (pre) SE (pre) M (unadj) SE (unadj) M (adj) SE (adj) conf.low conf.high
Controle F 208 40.635 1.278 38.615 1.351 35.106 0.503 34.119 36.093
Controle M 211 33.934 1.457 32.393 1.517 35.389 0.499 34.410 36.368
Experimental F 276 39.623 1.053 40.272 1.123 37.744 0.436 36.889 38.600
Experimental M 279 34.082 1.030 34.484 1.053 37.335 0.434 36.483 38.187

Plots for ancova

if (length(unique(pdat[["Sexo"]])) >= 2) {
  ggPlotAoC2(pwcs, "grupo", "Sexo", aov, ylab = "Writing (TDE)",
             subtitle = which(aov$Effect == "grupo:Sexo"), addParam = "errorbar") +
    ggplot2::scale_color_manual(values = color[["Sexo"]])  +
    ggplot2::ylab("Writing (TDE)") +
    theme(axis.title = element_text(size = 14),
          legend.text = element_text(size = 16),
          plot.subtitle = element_text(size = 18)) +
    if (ymin.ci < ymax.ci) ggplot2::ylim(ymin.ci, ymax.ci)
}
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.

if (length(unique(pdat[["Sexo"]])) >= 2) {
  ggPlotAoC2(pwcs, "Sexo", "grupo", aov, ylab = "Writing (TDE)",
               subtitle = which(aov$Effect == "grupo:Sexo"), addParam = "errorbar") +
      ggplot2::scale_color_manual(values = color[["grupo"]])  +
      ggplot2::ylab("Writing (TDE)") +
      theme(axis.title = element_text(size = 14),
            legend.text = element_text(size = 16),
            plot.subtitle = element_text(size = 18)) +
      if (ymin.ci < ymax.ci) ggplot2::ylim(ymin.ci, ymax.ci)
}
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.

if (length(unique(pdat[["Sexo"]])) >= 2) {
  plots <- twoWayAncovaBoxPlots(
    wdat, "score.tde.pos", c("grupo","Sexo"), aov, pwcs, covar = "score.tde.pre",
    theme = "classic", color = color[["grupo:Sexo"]],
    subtitle = which(aov$Effect == "grupo:Sexo"))
}
if (length(unique(pdat[["Sexo"]])) >= 2) {
  plots[["grupo:Sexo"]] + ggplot2::ylab("Writing (TDE)") +
  ggplot2::scale_x_discrete(labels=c('pre', 'pos')) +
  if (ymin < ymax) ggplot2::ylim(ymin, ymax)
}
## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's colour values.

if (length(unique(pdat[["Sexo"]])) >= 2) {
  plots <- twoWayAncovaBoxPlots(
    wdat.long, "score.tde", c("grupo","Sexo"), aov, pwc.long,
    pre.post = "time",
    theme = "classic", color = color$prepost)
}
if (length(unique(pdat[["Sexo"]])) >= 2) 
  plots[["grupo:Sexo"]] + ggplot2::ylab("Writing (TDE)") +
    if (ymin < ymax) ggplot2::ylim(ymin, ymax)

Checking linearity assumption

if (length(unique(pdat[["Sexo"]])) >= 2) {
  ggscatter(wdat, x = "score.tde.pre", y = "score.tde.pos", size = 0.5,
            facet.by = c("grupo","Sexo"), add = "reg.line")+
    stat_regline_equation(
      aes(label =  paste(..eq.label.., ..rr.label.., sep = "~~~~"))
    ) +
    ggplot2::xlab("Writing (TDE)") +
    ggplot2::ylab("Writing (TDE)") +
    theme(axis.title = element_text(size = 14),
          legend.text = element_text(size = 16),
          plot.subtitle = element_text(size = 18)) +
    if (ymin < ymax) ggplot2::ylim(ymin, ymax)
}

if (length(unique(pdat[["Sexo"]])) >= 2) {
  ggscatter(wdat, x = "score.tde.pre", y = "score.tde.pos", size = 0.5,
            color = "grupo", facet.by = "Sexo", add = "reg.line")+
    stat_regline_equation(
      aes(label =  paste(..eq.label.., ..rr.label.., sep = "~~~~"), color = grupo)
    ) +
    ggplot2::labs(subtitle = rstatix::get_test_label(aov, detailed = T, row = which(aov$Effect == "grupo:Sexo"))) +
    ggplot2::scale_color_manual(values = color[["grupo"]]) +
    ggplot2::xlab("Writing (TDE)") +
    ggplot2::ylab("Writing (TDE)") +
    theme(axis.title = element_text(size = 14),
          legend.text = element_text(size = 16),
          plot.subtitle = element_text(size = 18)) +
    if (ymin < ymax) ggplot2::ylim(ymin, ymax)
}

if (length(unique(pdat[["Sexo"]])) >= 2) {
  ggscatter(wdat, x = "score.tde.pre", y = "score.tde.pos", size = 0.5,
            color = "Sexo", facet.by = "grupo", add = "reg.line")+
    stat_regline_equation(
      aes(label =  paste(..eq.label.., ..rr.label.., sep = "~~~~"), color = Sexo)
    ) +
    ggplot2::labs(subtitle = rstatix::get_test_label(aov, detailed = T, row = which(aov$Effect == "grupo:Sexo"))) +
    ggplot2::scale_color_manual(values = color[["Sexo"]]) +
    ggplot2::xlab("Writing (TDE)") +
    ggplot2::ylab("Writing (TDE)") +
    theme(axis.title = element_text(size = 14),
          legend.text = element_text(size = 16),
          plot.subtitle = element_text(size = 18)) +
    if (ymin < ymax) ggplot2::ylim(ymin, ymax)
}

Checking normality and homogeneity

if (length(unique(pdat[["Sexo"]])) >= 2) 
  res <- augment(lm(score.tde.pos ~ score.tde.pre + grupo*Sexo, data = wdat))
if (length(unique(pdat[["Sexo"]])) >= 2)
  shapiro_test(res$.resid)
## # A tibble: 1 × 3
##   variable   statistic p.value
##   <chr>          <dbl>   <dbl>
## 1 res$.resid     0.997  0.0375
if (length(unique(pdat[["Sexo"]])) >= 2) 
  levene_test(res, .resid ~ grupo*Sexo)
## # A tibble: 1 × 4
##     df1   df2 statistic     p
##   <int> <int>     <dbl> <dbl>
## 1     3   970    0.0813 0.970

ANCOVA and Pairwise for two factors grupo:Zona

Without remove non-normal data

pdat = remove_group_data(dat[!is.na(dat[["grupo"]]) & !is.na(dat[["Zona"]]),],
                         "score.tde.pos", c("grupo","Zona"))
pdat = pdat[pdat[["Zona"]] %in% do.call(
  intersect, lapply(unique(pdat[["grupo"]]), FUN = function(x) {
    unique(pdat[["Zona"]][which(pdat[["grupo"]] == x)])
  })),]
pdat[["grupo"]] = factor(pdat[["grupo"]], level[["grupo"]])
pdat[["Zona"]] = factor(
  pdat[["Zona"]],
  level[["Zona"]][level[["Zona"]] %in% unique(pdat[["Zona"]])])

pdat.long <- rbind(pdat[,c("id","grupo","Zona")], pdat[,c("id","grupo","Zona")])
pdat.long[["time"]] <- c(rep("pre", nrow(pdat)), rep("pos", nrow(pdat)))
pdat.long[["time"]] <- factor(pdat.long[["time"]], c("pre","pos"))
pdat.long[["score.tde"]] <- c(pdat[["score.tde.pre"]], pdat[["score.tde.pos"]])

if (length(unique(pdat[["Zona"]])) >= 2) {
  aov = anova_test(pdat, score.tde.pos ~ score.tde.pre + grupo*Zona)
  laov[["grupo:Zona"]] <- get_anova_table(aov)
}
if (length(unique(pdat[["Zona"]])) >= 2) {
  pwcs <- list()
  pwcs[["Zona"]] <- emmeans_test(
    group_by(pdat, grupo), score.tde.pos ~ Zona,
    covariate = score.tde.pre, p.adjust.method = "bonferroni")
  pwcs[["grupo"]] <- emmeans_test(
    group_by(pdat, Zona), score.tde.pos ~ grupo,
    covariate = score.tde.pre, p.adjust.method = "bonferroni")
  
  pwc <- plyr::rbind.fill(pwcs[["grupo"]], pwcs[["Zona"]])
  pwc <- pwc[,c("grupo","Zona", colnames(pwc)[!colnames(pwc) %in% c("grupo","Zona")])]
}
if (length(unique(pdat[["Zona"]])) >= 2) {
  pwc.long <- emmeans_test(dplyr::group_by_at(pdat.long, c("grupo","Zona")),
                           score.tde ~ time,
                           p.adjust.method = "bonferroni")
  lpwc[["grupo:Zona"]] <- plyr::rbind.fill(pwc, pwc.long)
}
if (length(unique(pdat[["Zona"]])) >= 2) {
  ds <- get.descriptives(pdat, "score.tde.pos", c("grupo","Zona"), covar = "score.tde.pre")
  ds <- merge(ds[ds$variable != "score.tde.pre",],
              ds[ds$variable == "score.tde.pre", !colnames(ds) %in% c("variable")],
              by = c("grupo","Zona"), all.x = T, suffixes = c("", ".score.tde.pre"))
  ds <- merge(get_emmeans(pwcs[["grupo"]]), ds,
              by = c("grupo","Zona"), suffixes = c(".emms", ""))
  ds <- ds[,c("grupo","Zona","n","mean.score.tde.pre","se.score.tde.pre","mean","se",
              "emmean","se.emms","conf.low","conf.high")]
  
  colnames(ds) <- c("grupo","Zona", "N", paste0(c("M","SE")," (pre)"),
                    paste0(c("M","SE"), " (unadj)"),
                    paste0(c("M", "SE"), " (adj)"), "conf.low", "conf.high")
  
  lemms[["grupo:Zona"]] <- ds
}

Computing ANCOVA and PairWise After removing non-normal data (OK)

if (length(unique(pdat[["Zona"]])) >= 2) {
  wdat = pdat 
  
  res = residuals(lm(score.tde.pos ~ score.tde.pre + grupo*Zona, data = wdat))
  non.normal = getNonNormal(res, wdat$id, plimit = 0.05)
  
  wdat = wdat[!wdat$id %in% non.normal,]
  
  wdat.long <- rbind(wdat[,c("id","grupo","Zona")], wdat[,c("id","grupo","Zona")])
  wdat.long[["time"]] <- c(rep("pre", nrow(wdat)), rep("pos", nrow(wdat)))
  wdat.long[["time"]] <- factor(wdat.long[["time"]], c("pre","pos"))
  wdat.long[["score.tde"]] <- c(wdat[["score.tde.pre"]], wdat[["score.tde.pos"]])
  
  
  ldat[["grupo:Zona"]] = wdat
  
  (non.normal)
}
##  [1] "P962"  "P984"  "P1971" "P2929" "P1117" "P2848" "P2917" "P2870" "P2871"
## [10] "P2835" "P921"  "P2959" "P3029" "P1139" "P572"  "P2983" "P2858" "P3021"
## [19] "P3080" "P606"  "P2933" "P914"  "P2831" "P2946" "P3533" "P2953" "P2886"
## [28] "P3548" "P3476" "P2861" "P2880" "P2982" "P2995" "P3674" "P2864" "P2883"
## [37] "P2950" "P2004" "P3475" "P3012" "P3010" "P3015" "P2975" "P2937" "P3660"
## [46] "P2910" "P1691" "P2868" "P2860" "P1878" "P609"  "P908"  "P2997" "P2956"
## [55] "P2973" "P2879" "P2854" "P3024" "P2967" "P2948" "P2909" "P1056" "P2964"
if (length(unique(pdat[["Zona"]])) >= 2) {
  aov = anova_test(wdat, score.tde.pos ~ score.tde.pre + grupo*Zona)
  laov[["grupo:Zona"]] <- merge(get_anova_table(aov), laov[["grupo:Zona"]],
                                         by="Effect", suffixes = c("","'"))
  df = get_anova_table(aov)
}
Effect DFn DFd F p p<.05 ges
score.tde.pre 1 735 3559.226 0 * 0.829
grupo 1 735 26.138 0 * 0.034
Zona 1 735 54.679 0 * 0.069
grupo:Zona 1 735 14.155 0 * 0.019 
if (length(unique(pdat[["Zona"]])) >= 2) {
  pwcs <- list()
  pwcs[["Zona"]] <- emmeans_test(
    group_by(wdat, grupo), score.tde.pos ~ Zona,
    covariate = score.tde.pre, p.adjust.method = "bonferroni")
  pwcs[["grupo"]] <- emmeans_test(
    group_by(wdat, Zona), score.tde.pos ~ grupo,
    covariate = score.tde.pre, p.adjust.method = "bonferroni")
  
  pwc <- plyr::rbind.fill(pwcs[["grupo"]], pwcs[["Zona"]])
  pwc <- pwc[,c("grupo","Zona", colnames(pwc)[!colnames(pwc) %in% c("grupo","Zona")])]
}
grupo Zona term .y. group1 group2 df statistic p p.adj p.adj.signif
Rural score.tde.pre*grupo score.tde.pos Controle Experimental 735 -2.009 0.045 0.045 *
Urbana score.tde.pre*grupo score.tde.pos Controle Experimental 735 -6.023 0.000 0.000 ****
Controle score.tde.pre*Zona score.tde.pos Rural Urbana 735 7.675 0.000 0.000 ****
Experimental score.tde.pre*Zona score.tde.pos Rural Urbana 735 3.188 0.001 0.001 ** 
if (length(unique(pdat[["Zona"]])) >= 2) {
  pwc.long <- emmeans_test(dplyr::group_by_at(wdat.long, c("grupo","Zona")),
                           score.tde ~ time,
                           p.adjust.method = "bonferroni")
  lpwc[["grupo:Zona"]] <- merge(plyr::rbind.fill(pwc, pwc.long),
                                         lpwc[["grupo:Zona"]],
                                         by=c("grupo","Zona","term",".y.","group1","group2"),
                                         suffixes = c("","'"))
}
grupo Zona term .y. group1 group2 df statistic p p.adj p.adj.signif
Controle Rural time score.tde pre pos 1472 0.543 0.588 0.588 ns
Controle Urbana time score.tde pre pos 1472 3.191 0.001 0.001 **
Experimental Rural time score.tde pre pos 1472 -0.333 0.739 0.739 ns
Experimental Urbana time score.tde pre pos 1472 1.001 0.317 0.317 ns 
if (length(unique(pdat[["Zona"]])) >= 2) {
  ds <- get.descriptives(wdat, "score.tde.pos", c("grupo","Zona"), covar = "score.tde.pre")
  ds <- merge(ds[ds$variable != "score.tde.pre",],
              ds[ds$variable == "score.tde.pre", !colnames(ds) %in% c("variable")],
              by = c("grupo","Zona"), all.x = T, suffixes = c("", ".score.tde.pre"))
  ds <- merge(get_emmeans(pwcs[["grupo"]]), ds,
              by = c("grupo","Zona"), suffixes = c(".emms", ""))
  ds <- ds[,c("grupo","Zona","n","mean.score.tde.pre","se.score.tde.pre",
              "mean","se","emmean","se.emms","conf.low","conf.high")]
  
  colnames(ds) <- c("grupo","Zona", "N", paste0(c("M","SE")," (pre)"),
                    paste0(c("M","SE"), " (unadj)"),
                    paste0(c("M", "SE"), " (adj)"), "conf.low", "conf.high")
  
  lemms[["grupo:Zona"]] <- merge(ds, lemms[["grupo:Zona"]],
                                          by=c("grupo","Zona"), suffixes = c("","'"))
}
grupo Zona N M (pre) SE (pre) M (unadj) SE (unadj) M (adj) SE (adj) conf.low conf.high
Controle Rural 228 35.956 1.330 34.978 1.426 35.358 0.547 34.284 36.432
Controle Urbana 101 39.743 1.823 31.099 2.044 27.763 0.824 26.145 29.380
Experimental Rural 260 34.215 1.113 34.777 1.186 36.865 0.514 35.857 37.874
Experimental Urbana 151 38.318 1.389 36.099 1.501 34.161 0.673 32.840 35.483

Plots for ancova

if (length(unique(pdat[["Zona"]])) >= 2) {
  ggPlotAoC2(pwcs, "grupo", "Zona", aov, ylab = "Writing (TDE)",
             subtitle = which(aov$Effect == "grupo:Zona"), addParam = "errorbar") +
    ggplot2::scale_color_manual(values = color[["Zona"]])  +
    ggplot2::ylab("Writing (TDE)") +
    theme(axis.title = element_text(size = 14),
          legend.text = element_text(size = 16),
          plot.subtitle = element_text(size = 18)) +
    if (ymin.ci < ymax.ci) ggplot2::ylim(ymin.ci, ymax.ci)
}
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.

if (length(unique(pdat[["Zona"]])) >= 2) {
  ggPlotAoC2(pwcs, "Zona", "grupo", aov, ylab = "Writing (TDE)",
               subtitle = which(aov$Effect == "grupo:Zona"), addParam = "errorbar") +
      ggplot2::scale_color_manual(values = color[["grupo"]])  +
      ggplot2::ylab("Writing (TDE)") +
      theme(axis.title = element_text(size = 14),
            legend.text = element_text(size = 16),
            plot.subtitle = element_text(size = 18)) +
      if (ymin.ci < ymax.ci) ggplot2::ylim(ymin.ci, ymax.ci)
}
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.

if (length(unique(pdat[["Zona"]])) >= 2) {
  plots <- twoWayAncovaBoxPlots(
    wdat, "score.tde.pos", c("grupo","Zona"), aov, pwcs, covar = "score.tde.pre",
    theme = "classic", color = color[["grupo:Zona"]],
    subtitle = which(aov$Effect == "grupo:Zona"))
}
if (length(unique(pdat[["Zona"]])) >= 2) {
  plots[["grupo:Zona"]] + ggplot2::ylab("Writing (TDE)") +
  ggplot2::scale_x_discrete(labels=c('pre', 'pos')) +
  if (ymin < ymax) ggplot2::ylim(ymin, ymax)
}
## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's colour values.

if (length(unique(pdat[["Zona"]])) >= 2) {
  plots <- twoWayAncovaBoxPlots(
    wdat.long, "score.tde", c("grupo","Zona"), aov, pwc.long,
    pre.post = "time",
    theme = "classic", color = color$prepost)
}
if (length(unique(pdat[["Zona"]])) >= 2) 
  plots[["grupo:Zona"]] + ggplot2::ylab("Writing (TDE)") +
    if (ymin < ymax) ggplot2::ylim(ymin, ymax)

Checking linearity assumption

if (length(unique(pdat[["Zona"]])) >= 2) {
  ggscatter(wdat, x = "score.tde.pre", y = "score.tde.pos", size = 0.5,
            facet.by = c("grupo","Zona"), add = "reg.line")+
    stat_regline_equation(
      aes(label =  paste(..eq.label.., ..rr.label.., sep = "~~~~"))
    ) +
    ggplot2::xlab("Writing (TDE)") +
    ggplot2::ylab("Writing (TDE)") +
    theme(axis.title = element_text(size = 14),
          legend.text = element_text(size = 16),
          plot.subtitle = element_text(size = 18)) +
    if (ymin < ymax) ggplot2::ylim(ymin, ymax)
}

if (length(unique(pdat[["Zona"]])) >= 2) {
  ggscatter(wdat, x = "score.tde.pre", y = "score.tde.pos", size = 0.5,
            color = "grupo", facet.by = "Zona", add = "reg.line")+
    stat_regline_equation(
      aes(label =  paste(..eq.label.., ..rr.label.., sep = "~~~~"), color = grupo)
    ) +
    ggplot2::labs(subtitle = rstatix::get_test_label(aov, detailed = T, row = which(aov$Effect == "grupo:Zona"))) +
    ggplot2::scale_color_manual(values = color[["grupo"]]) +
    ggplot2::xlab("Writing (TDE)") +
    ggplot2::ylab("Writing (TDE)") +
    theme(axis.title = element_text(size = 14),
          legend.text = element_text(size = 16),
          plot.subtitle = element_text(size = 18)) +
    if (ymin < ymax) ggplot2::ylim(ymin, ymax)
}

if (length(unique(pdat[["Zona"]])) >= 2) {
  ggscatter(wdat, x = "score.tde.pre", y = "score.tde.pos", size = 0.5,
            color = "Zona", facet.by = "grupo", add = "reg.line")+
    stat_regline_equation(
      aes(label =  paste(..eq.label.., ..rr.label.., sep = "~~~~"), color = Zona)
    ) +
    ggplot2::labs(subtitle = rstatix::get_test_label(aov, detailed = T, row = which(aov$Effect == "grupo:Zona"))) +
    ggplot2::scale_color_manual(values = color[["Zona"]]) +
    ggplot2::xlab("Writing (TDE)") +
    ggplot2::ylab("Writing (TDE)") +
    theme(axis.title = element_text(size = 14),
          legend.text = element_text(size = 16),
          plot.subtitle = element_text(size = 18)) +
    if (ymin < ymax) ggplot2::ylim(ymin, ymax)
}

Checking normality and homogeneity

if (length(unique(pdat[["Zona"]])) >= 2) 
  res <- augment(lm(score.tde.pos ~ score.tde.pre + grupo*Zona, data = wdat))
if (length(unique(pdat[["Zona"]])) >= 2)
  shapiro_test(res$.resid)
## # A tibble: 1 × 3
##   variable   statistic p.value
##   <chr>          <dbl>   <dbl>
## 1 res$.resid     0.995  0.0253
if (length(unique(pdat[["Zona"]])) >= 2) 
  levene_test(res, .resid ~ grupo*Zona)
## # A tibble: 1 × 4
##     df1   df2 statistic      p
##   <int> <int>     <dbl>  <dbl>
## 1     3   736      2.40 0.0667

ANCOVA and Pairwise for two factors grupo:Cor.Raca

Without remove non-normal data

pdat = remove_group_data(dat[!is.na(dat[["grupo"]]) & !is.na(dat[["Cor.Raca"]]),],
                         "score.tde.pos", c("grupo","Cor.Raca"))
## Warning: There was 1 warning in `mutate()`.
## ℹ In argument: `ci = abs(stats::qt(alpha/2, .data$n - 1) * .data$se)`.
## Caused by warning:
## ! There was 1 warning in `mutate()`.
## ℹ In argument: `ci = abs(stats::qt(alpha/2, .data$n - 1) * .data$se)`.
## Caused by warning in `stats::qt()`:
## ! NaNs produced
pdat = pdat[pdat[["Cor.Raca"]] %in% do.call(
  intersect, lapply(unique(pdat[["grupo"]]), FUN = function(x) {
    unique(pdat[["Cor.Raca"]][which(pdat[["grupo"]] == x)])
  })),]
pdat[["grupo"]] = factor(pdat[["grupo"]], level[["grupo"]])
pdat[["Cor.Raca"]] = factor(
  pdat[["Cor.Raca"]],
  level[["Cor.Raca"]][level[["Cor.Raca"]] %in% unique(pdat[["Cor.Raca"]])])

pdat.long <- rbind(pdat[,c("id","grupo","Cor.Raca")], pdat[,c("id","grupo","Cor.Raca")])
pdat.long[["time"]] <- c(rep("pre", nrow(pdat)), rep("pos", nrow(pdat)))
pdat.long[["time"]] <- factor(pdat.long[["time"]], c("pre","pos"))
pdat.long[["score.tde"]] <- c(pdat[["score.tde.pre"]], pdat[["score.tde.pos"]])

if (length(unique(pdat[["Cor.Raca"]])) >= 2) {
  aov = anova_test(pdat, score.tde.pos ~ score.tde.pre + grupo*Cor.Raca)
  laov[["grupo:Cor.Raca"]] <- get_anova_table(aov)
}
if (length(unique(pdat[["Cor.Raca"]])) >= 2) {
  pwcs <- list()
  pwcs[["Cor.Raca"]] <- emmeans_test(
    group_by(pdat, grupo), score.tde.pos ~ Cor.Raca,
    covariate = score.tde.pre, p.adjust.method = "bonferroni")
  pwcs[["grupo"]] <- emmeans_test(
    group_by(pdat, Cor.Raca), score.tde.pos ~ grupo,
    covariate = score.tde.pre, p.adjust.method = "bonferroni")
  
  pwc <- plyr::rbind.fill(pwcs[["grupo"]], pwcs[["Cor.Raca"]])
  pwc <- pwc[,c("grupo","Cor.Raca", colnames(pwc)[!colnames(pwc) %in% c("grupo","Cor.Raca")])]
}
if (length(unique(pdat[["Cor.Raca"]])) >= 2) {
  pwc.long <- emmeans_test(dplyr::group_by_at(pdat.long, c("grupo","Cor.Raca")),
                           score.tde ~ time,
                           p.adjust.method = "bonferroni")
  lpwc[["grupo:Cor.Raca"]] <- plyr::rbind.fill(pwc, pwc.long)
}
if (length(unique(pdat[["Cor.Raca"]])) >= 2) {
  ds <- get.descriptives(pdat, "score.tde.pos", c("grupo","Cor.Raca"), covar = "score.tde.pre")
  ds <- merge(ds[ds$variable != "score.tde.pre",],
              ds[ds$variable == "score.tde.pre", !colnames(ds) %in% c("variable")],
              by = c("grupo","Cor.Raca"), all.x = T, suffixes = c("", ".score.tde.pre"))
  ds <- merge(get_emmeans(pwcs[["grupo"]]), ds,
              by = c("grupo","Cor.Raca"), suffixes = c(".emms", ""))
  ds <- ds[,c("grupo","Cor.Raca","n","mean.score.tde.pre","se.score.tde.pre","mean","se",
              "emmean","se.emms","conf.low","conf.high")]
  
  colnames(ds) <- c("grupo","Cor.Raca", "N", paste0(c("M","SE")," (pre)"),
                    paste0(c("M","SE"), " (unadj)"),
                    paste0(c("M", "SE"), " (adj)"), "conf.low", "conf.high")
  
  lemms[["grupo:Cor.Raca"]] <- ds
}

Computing ANCOVA and PairWise After removing non-normal data (OK)

if (length(unique(pdat[["Cor.Raca"]])) >= 2) {
  wdat = pdat 
  
  res = residuals(lm(score.tde.pos ~ score.tde.pre + grupo*Cor.Raca, data = wdat))
  non.normal = getNonNormal(res, wdat$id, plimit = 0.05)
  
  wdat = wdat[!wdat$id %in% non.normal,]
  
  wdat.long <- rbind(wdat[,c("id","grupo","Cor.Raca")], wdat[,c("id","grupo","Cor.Raca")])
  wdat.long[["time"]] <- c(rep("pre", nrow(wdat)), rep("pos", nrow(wdat)))
  wdat.long[["time"]] <- factor(wdat.long[["time"]], c("pre","pos"))
  wdat.long[["score.tde"]] <- c(wdat[["score.tde.pre"]], wdat[["score.tde.pos"]])
  
  
  ldat[["grupo:Cor.Raca"]] = wdat
  
  (non.normal)
}
##  [1] "P1128" "P1126" "P908"  "P2995" "P1139" "P3021" "P2871" "P2858" "P3637"
## [10] "P2870" "P1117" "P572"  "P2978" "P2975" "P2983" "P2904" "P2835" "P3000"
## [19] "P2888" "P2993" "P2905" "P1118" "P1804" "P3007" "P3228" "P3029" "P2865"
## [28] "P2967" "P1018" "P2994" "P2953" "P2997" "P3476" "P3054" "P606"
if (length(unique(pdat[["Cor.Raca"]])) >= 2) {
  aov = anova_test(wdat, score.tde.pos ~ score.tde.pre + grupo*Cor.Raca)
  laov[["grupo:Cor.Raca"]] <- merge(get_anova_table(aov), laov[["grupo:Cor.Raca"]],
                                         by="Effect", suffixes = c("","'"))
  df = get_anova_table(aov)
}
Effect DFn DFd F p p<.05 ges
score.tde.pre 1 447 1620.497 0.000 * 0.784
grupo 1 447 7.311 0.007 * 0.016
Cor.Raca 2 447 3.357 0.036 * 0.015
grupo:Cor.Raca 2 447 1.174 0.310 0.005 
if (length(unique(pdat[["Cor.Raca"]])) >= 2) {
  pwcs <- list()
  pwcs[["Cor.Raca"]] <- emmeans_test(
    group_by(wdat, grupo), score.tde.pos ~ Cor.Raca,
    covariate = score.tde.pre, p.adjust.method = "bonferroni")
  pwcs[["grupo"]] <- emmeans_test(
    group_by(wdat, Cor.Raca), score.tde.pos ~ grupo,
    covariate = score.tde.pre, p.adjust.method = "bonferroni")
  
  pwc <- plyr::rbind.fill(pwcs[["grupo"]], pwcs[["Cor.Raca"]])
  pwc <- pwc[,c("grupo","Cor.Raca", colnames(pwc)[!colnames(pwc) %in% c("grupo","Cor.Raca")])]
}
grupo Cor.Raca term .y. group1 group2 df statistic p p.adj p.adj.signif
Parda score.tde.pre*grupo score.tde.pos Controle Experimental 447 -1.911 0.057 0.057 ns
Indígena score.tde.pre*grupo score.tde.pos Controle Experimental 447 0.274 0.784 0.784 ns
Branca score.tde.pre*grupo score.tde.pos Controle Experimental 447 -2.442 0.015 0.015 *
Controle score.tde.pre*Cor.Raca score.tde.pos Parda Indígena 447 -2.163 0.031 0.093 ns
Controle score.tde.pre*Cor.Raca score.tde.pos Parda Branca 447 1.334 0.183 0.549 ns
Controle score.tde.pre*Cor.Raca score.tde.pos Indígena Branca 447 2.687 0.007 0.022 *
Experimental score.tde.pre*Cor.Raca score.tde.pos Parda Indígena 447 -1.314 0.190 0.569 ns
Experimental score.tde.pre*Cor.Raca score.tde.pos Parda Branca 447 -0.308 0.759 1.000 ns
Experimental score.tde.pre*Cor.Raca score.tde.pos Indígena Branca 447 1.060 0.290 0.870 ns 
if (length(unique(pdat[["Cor.Raca"]])) >= 2) {
  pwc.long <- emmeans_test(dplyr::group_by_at(wdat.long, c("grupo","Cor.Raca")),
                           score.tde ~ time,
                           p.adjust.method = "bonferroni")
  lpwc[["grupo:Cor.Raca"]] <- merge(plyr::rbind.fill(pwc, pwc.long),
                                         lpwc[["grupo:Cor.Raca"]],
                                         by=c("grupo","Cor.Raca","term",".y.","group1","group2"),
                                         suffixes = c("","'"))
}
grupo Cor.Raca term .y. group1 group2 df statistic p p.adj p.adj.signif
Controle Parda time score.tde pre pos 896 1.090 0.276 0.276 ns
Controle Indígena time score.tde pre pos 896 -0.441 0.659 0.659 ns
Controle Branca time score.tde pre pos 896 1.172 0.241 0.241 ns
Experimental Parda time score.tde pre pos 896 0.163 0.870 0.870 ns
Experimental Indígena time score.tde pre pos 896 -0.463 0.644 0.644 ns
Experimental Branca time score.tde pre pos 896 -0.019 0.985 0.985 ns 
if (length(unique(pdat[["Cor.Raca"]])) >= 2) {
  ds <- get.descriptives(wdat, "score.tde.pos", c("grupo","Cor.Raca"), covar = "score.tde.pre")
  ds <- merge(ds[ds$variable != "score.tde.pre",],
              ds[ds$variable == "score.tde.pre", !colnames(ds) %in% c("variable")],
              by = c("grupo","Cor.Raca"), all.x = T, suffixes = c("", ".score.tde.pre"))
  ds <- merge(get_emmeans(pwcs[["grupo"]]), ds,
              by = c("grupo","Cor.Raca"), suffixes = c(".emms", ""))
  ds <- ds[,c("grupo","Cor.Raca","n","mean.score.tde.pre","se.score.tde.pre",
              "mean","se","emmean","se.emms","conf.low","conf.high")]
  
  colnames(ds) <- c("grupo","Cor.Raca", "N", paste0(c("M","SE")," (pre)"),
                    paste0(c("M","SE"), " (unadj)"),
                    paste0(c("M", "SE"), " (adj)"), "conf.low", "conf.high")
  
  lemms[["grupo:Cor.Raca"]] <- merge(ds, lemms[["grupo:Cor.Raca"]],
                                          by=c("grupo","Cor.Raca"), suffixes = c("","'"))
}
grupo Cor.Raca N M (pre) SE (pre) M (unadj) SE (unadj) M (adj) SE (adj) conf.low conf.high
Controle Branca 45 41.378 2.626 36.600 3.007 30.951 1.399 28.202 33.701
Controle Indígena 11 42.000 5.215 45.636 6.009 39.393 2.820 33.851 44.935
Controle Parda 150 36.320 1.578 33.887 1.695 33.075 0.763 31.576 34.574
Experimental Branca 58 34.517 2.651 34.586 2.746 35.499 1.226 33.088 37.909
Experimental Indígena 15 27.067 4.406 30.333 4.574 38.371 2.420 33.616 43.126
Experimental Parda 175 33.851 1.343 33.514 1.460 35.064 0.707 33.674 36.453

Plots for ancova

if (length(unique(pdat[["Cor.Raca"]])) >= 2) {
  ggPlotAoC2(pwcs, "grupo", "Cor.Raca", aov, ylab = "Writing (TDE)",
             subtitle = which(aov$Effect == "grupo:Cor.Raca"), addParam = "errorbar") +
    ggplot2::scale_color_manual(values = color[["Cor.Raca"]])  +
    ggplot2::ylab("Writing (TDE)") +
    theme(axis.title = element_text(size = 14),
          legend.text = element_text(size = 16),
          plot.subtitle = element_text(size = 18)) +
    if (ymin.ci < ymax.ci) ggplot2::ylim(ymin.ci, ymax.ci)
}
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.

if (length(unique(pdat[["Cor.Raca"]])) >= 2) {
  ggPlotAoC2(pwcs, "Cor.Raca", "grupo", aov, ylab = "Writing (TDE)",
               subtitle = which(aov$Effect == "grupo:Cor.Raca"), addParam = "errorbar") +
      ggplot2::scale_color_manual(values = color[["grupo"]])  +
      ggplot2::ylab("Writing (TDE)") +
      theme(axis.title = element_text(size = 14),
            legend.text = element_text(size = 16),
            plot.subtitle = element_text(size = 18)) +
      if (ymin.ci < ymax.ci) ggplot2::ylim(ymin.ci, ymax.ci)
}
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.

if (length(unique(pdat[["Cor.Raca"]])) >= 2) {
  plots <- twoWayAncovaBoxPlots(
    wdat, "score.tde.pos", c("grupo","Cor.Raca"), aov, pwcs, covar = "score.tde.pre",
    theme = "classic", color = color[["grupo:Cor.Raca"]],
    subtitle = which(aov$Effect == "grupo:Cor.Raca"))
}
if (length(unique(pdat[["Cor.Raca"]])) >= 2) {
  plots[["grupo:Cor.Raca"]] + ggplot2::ylab("Writing (TDE)") +
  ggplot2::scale_x_discrete(labels=c('pre', 'pos')) +
  if (ymin < ymax) ggplot2::ylim(ymin, ymax)
}
## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's colour values.

if (length(unique(pdat[["Cor.Raca"]])) >= 2) {
  plots <- twoWayAncovaBoxPlots(
    wdat.long, "score.tde", c("grupo","Cor.Raca"), aov, pwc.long,
    pre.post = "time",
    theme = "classic", color = color$prepost)
}
if (length(unique(pdat[["Cor.Raca"]])) >= 2) 
  plots[["grupo:Cor.Raca"]] + ggplot2::ylab("Writing (TDE)") +
    if (ymin < ymax) ggplot2::ylim(ymin, ymax)

Checking linearity assumption

if (length(unique(pdat[["Cor.Raca"]])) >= 2) {
  ggscatter(wdat, x = "score.tde.pre", y = "score.tde.pos", size = 0.5,
            facet.by = c("grupo","Cor.Raca"), add = "reg.line")+
    stat_regline_equation(
      aes(label =  paste(..eq.label.., ..rr.label.., sep = "~~~~"))
    ) +
    ggplot2::xlab("Writing (TDE)") +
    ggplot2::ylab("Writing (TDE)") +
    theme(axis.title = element_text(size = 14),
        legend.text = element_text(size = 16),
        plot.subtitle = element_text(size = 18)) +
    if (ymin < ymax) ggplot2::ylim(ymin, ymax)
}

if (length(unique(pdat[["Cor.Raca"]])) >= 2) {
  ggscatter(wdat, x = "score.tde.pre", y = "score.tde.pos", size = 0.5,
            color = "grupo", facet.by = "Cor.Raca", add = "reg.line")+
    stat_regline_equation(
      aes(label =  paste(..eq.label.., ..rr.label.., sep = "~~~~"), color = grupo)
    ) +
    ggplot2::labs(subtitle = rstatix::get_test_label(aov, detailed = T, row = which(aov$Effect == "grupo:Cor.Raca"))) +
    ggplot2::scale_color_manual(values = color[["grupo"]]) +
    ggplot2::xlab("Writing (TDE)") +
    ggplot2::ylab("Writing (TDE)") +
    theme(axis.title = element_text(size = 14),
          legend.text = element_text(size = 16),
          plot.subtitle = element_text(size = 18)) +
    if (ymin < ymax) ggplot2::ylim(ymin, ymax)
}

if (length(unique(pdat[["Cor.Raca"]])) >= 2) {
  ggscatter(wdat, x = "score.tde.pre", y = "score.tde.pos", size = 0.5,
            color = "Cor.Raca", facet.by = "grupo", add = "reg.line")+
    stat_regline_equation(
      aes(label =  paste(..eq.label.., ..rr.label.., sep = "~~~~"), color = Cor.Raca)
    ) +
    ggplot2::labs(subtitle = rstatix::get_test_label(aov, detailed = T, row = which(aov$Effect == "grupo:Cor.Raca"))) +
    ggplot2::scale_color_manual(values = color[["Cor.Raca"]]) +
    ggplot2::xlab("Writing (TDE)") +
    ggplot2::ylab("Writing (TDE)") +
    theme(axis.title = element_text(size = 14),
          legend.text = element_text(size = 16),
          plot.subtitle = element_text(size = 18)) +
    if (ymin < ymax) ggplot2::ylim(ymin, ymax)
}

Checking normality and homogeneity

if (length(unique(pdat[["Cor.Raca"]])) >= 2) 
  res <- augment(lm(score.tde.pos ~ score.tde.pre + grupo*Cor.Raca, data = wdat))
if (length(unique(pdat[["Cor.Raca"]])) >= 2)
  shapiro_test(res$.resid)
## # A tibble: 1 × 3
##   variable   statistic p.value
##   <chr>          <dbl>   <dbl>
## 1 res$.resid     0.992  0.0140
if (length(unique(pdat[["Cor.Raca"]])) >= 2) 
  levene_test(res, .resid ~ grupo*Cor.Raca)
## # A tibble: 1 × 4
##     df1   df2 statistic     p
##   <int> <int>     <dbl> <dbl>
## 1     5   448     0.202 0.962

ANCOVA and Pairwise for two factors grupo:Serie

Without remove non-normal data

pdat = remove_group_data(dat[!is.na(dat[["grupo"]]) & !is.na(dat[["Serie"]]),],
                         "score.tde.pos", c("grupo","Serie"))
pdat = pdat[pdat[["Serie"]] %in% do.call(
  intersect, lapply(unique(pdat[["grupo"]]), FUN = function(x) {
    unique(pdat[["Serie"]][which(pdat[["grupo"]] == x)])
  })),]
pdat[["grupo"]] = factor(pdat[["grupo"]], level[["grupo"]])
pdat[["Serie"]] = factor(
  pdat[["Serie"]],
  level[["Serie"]][level[["Serie"]] %in% unique(pdat[["Serie"]])])

pdat.long <- rbind(pdat[,c("id","grupo","Serie")], pdat[,c("id","grupo","Serie")])
pdat.long[["time"]] <- c(rep("pre", nrow(pdat)), rep("pos", nrow(pdat)))
pdat.long[["time"]] <- factor(pdat.long[["time"]], c("pre","pos"))
pdat.long[["score.tde"]] <- c(pdat[["score.tde.pre"]], pdat[["score.tde.pos"]])

if (length(unique(pdat[["Serie"]])) >= 2) {
  aov = anova_test(pdat, score.tde.pos ~ score.tde.pre + grupo*Serie)
  laov[["grupo:Serie"]] <- get_anova_table(aov)
}
if (length(unique(pdat[["Serie"]])) >= 2) {
  pwcs <- list()
  pwcs[["Serie"]] <- emmeans_test(
    group_by(pdat, grupo), score.tde.pos ~ Serie,
    covariate = score.tde.pre, p.adjust.method = "bonferroni")
  pwcs[["grupo"]] <- emmeans_test(
    group_by(pdat, Serie), score.tde.pos ~ grupo,
    covariate = score.tde.pre, p.adjust.method = "bonferroni")
  
  pwc <- plyr::rbind.fill(pwcs[["grupo"]], pwcs[["Serie"]])
  pwc <- pwc[,c("grupo","Serie", colnames(pwc)[!colnames(pwc) %in% c("grupo","Serie")])]
}
if (length(unique(pdat[["Serie"]])) >= 2) {
  pwc.long <- emmeans_test(dplyr::group_by_at(pdat.long, c("grupo","Serie")),
                           score.tde ~ time,
                           p.adjust.method = "bonferroni")
  lpwc[["grupo:Serie"]] <- plyr::rbind.fill(pwc, pwc.long)
}
if (length(unique(pdat[["Serie"]])) >= 2) {
  ds <- get.descriptives(pdat, "score.tde.pos", c("grupo","Serie"), covar = "score.tde.pre")
  ds <- merge(ds[ds$variable != "score.tde.pre",],
              ds[ds$variable == "score.tde.pre", !colnames(ds) %in% c("variable")],
              by = c("grupo","Serie"), all.x = T, suffixes = c("", ".score.tde.pre"))
  ds <- merge(get_emmeans(pwcs[["grupo"]]), ds,
              by = c("grupo","Serie"), suffixes = c(".emms", ""))
  ds <- ds[,c("grupo","Serie","n","mean.score.tde.pre","se.score.tde.pre","mean","se",
              "emmean","se.emms","conf.low","conf.high")]
  
  colnames(ds) <- c("grupo","Serie", "N", paste0(c("M","SE")," (pre)"),
                    paste0(c("M","SE"), " (unadj)"),
                    paste0(c("M", "SE"), " (adj)"), "conf.low", "conf.high")
  
  lemms[["grupo:Serie"]] <- ds
}

Computing ANCOVA and PairWise After removing non-normal data (OK)

if (length(unique(pdat[["Serie"]])) >= 2) {
  wdat = pdat 
  
  res = residuals(lm(score.tde.pos ~ score.tde.pre + grupo*Serie, data = wdat))
  non.normal = getNonNormal(res, wdat$id, plimit = 0.05)
  
  wdat = wdat[!wdat$id %in% non.normal,]
  
  wdat.long <- rbind(wdat[,c("id","grupo","Serie")], wdat[,c("id","grupo","Serie")])
  wdat.long[["time"]] <- c(rep("pre", nrow(wdat)), rep("pos", nrow(wdat)))
  wdat.long[["time"]] <- factor(wdat.long[["time"]], c("pre","pos"))
  wdat.long[["score.tde"]] <- c(wdat[["score.tde.pre"]], wdat[["score.tde.pos"]])
  
  
  ldat[["grupo:Serie"]] = wdat
  
  (non.normal)
}
##  [1] "P962"  "P1128" "P984"  "P1129" "P1971" "P1126" "P1139" "P2835" "P2983"
## [10] "P572"  "P1117" "P3571" "P3537" "P3574" "P1983" "P2848" "P809"  "P3637"
## [19] "P2910" "P976"  "P2870" "P2969" "P2866" "P3732" "P2913" "P921"  "P1018"
## [28] "P2964" "P1111" "P2967" "P3228" "P929"  "P2959" "P2886" "P914"  "P2995"
## [37] "P2982" "P2375" "P3241" "P2975" "P908"  "P2905" "P2986" "P1068" "P2946"
## [46] "P2904" "P1878" "P2883" "P2973" "P3533" "P3019" "P2865" "P2864" "P1056"
## [55] "P1152" "P2937" "P2876" "P3007" "P2978" "P2950" "P2871" "P2953" "P2917"
## [64] "P2948" "P609"
if (length(unique(pdat[["Serie"]])) >= 2) {
  aov = anova_test(wdat, score.tde.pos ~ score.tde.pre + grupo*Serie)
  laov[["grupo:Serie"]] <- merge(get_anova_table(aov), laov[["grupo:Serie"]],
                                         by="Effect", suffixes = c("","'"))
  df = get_anova_table(aov)
}
Effect DFn DFd F p p<.05 ges
score.tde.pre 1 1047 3720.788 0.000 * 0.780
grupo 1 1047 23.855 0.000 * 0.022
Serie 3 1047 34.727 0.000 * 0.090
grupo:Serie 3 1047 3.528 0.015 * 0.010 
if (length(unique(pdat[["Serie"]])) >= 2) {
  pwcs <- list()
  pwcs[["Serie"]] <- emmeans_test(
    group_by(wdat, grupo), score.tde.pos ~ Serie,
    covariate = score.tde.pre, p.adjust.method = "bonferroni")
  pwcs[["grupo"]] <- emmeans_test(
    group_by(wdat, Serie), score.tde.pos ~ grupo,
    covariate = score.tde.pre, p.adjust.method = "bonferroni")
  
  pwc <- plyr::rbind.fill(pwcs[["grupo"]], pwcs[["Serie"]])
  pwc <- pwc[,c("grupo","Serie", colnames(pwc)[!colnames(pwc) %in% c("grupo","Serie")])]
}
grupo Serie term .y. group1 group2 df statistic p p.adj p.adj.signif
6 ano score.tde.pre*grupo score.tde.pos Controle Experimental 1047 -0.249 0.803 0.803 ns
7 ano score.tde.pre*grupo score.tde.pos Controle Experimental 1047 -4.500 0.000 0.000 ****
8 ano score.tde.pre*grupo score.tde.pos Controle Experimental 1047 -3.348 0.001 0.001 ***
9 ano score.tde.pre*grupo score.tde.pos Controle Experimental 1047 -1.694 0.091 0.091 ns
Controle score.tde.pre*Serie score.tde.pos 6 ano 7 ano 1047 1.644 0.100 0.602 ns
Controle score.tde.pre*Serie score.tde.pos 6 ano 8 ano 1047 -2.433 0.015 0.091 ns
Controle score.tde.pre*Serie score.tde.pos 6 ano 9 ano 1047 -5.074 0.000 0.000 ****
Controle score.tde.pre*Serie score.tde.pos 7 ano 8 ano 1047 -3.939 0.000 0.001 ***
Controle score.tde.pre*Serie score.tde.pos 7 ano 9 ano 1047 -6.756 0.000 0.000 ****
Controle score.tde.pre*Serie score.tde.pos 8 ano 9 ano 1047 -2.300 0.022 0.130 ns
Experimental score.tde.pre*Serie score.tde.pos 6 ano 7 ano 1047 -2.574 0.010 0.061 ns
Experimental score.tde.pre*Serie score.tde.pos 6 ano 8 ano 1047 -6.514 0.000 0.000 ****
Experimental score.tde.pre*Serie score.tde.pos 6 ano 9 ano 1047 -7.114 0.000 0.000 ****
Experimental score.tde.pre*Serie score.tde.pos 7 ano 8 ano 1047 -4.242 0.000 0.000 ***
Experimental score.tde.pre*Serie score.tde.pos 7 ano 9 ano 1047 -4.930 0.000 0.000 ****
Experimental score.tde.pre*Serie score.tde.pos 8 ano 9 ano 1047 -0.688 0.492 1.000 ns 
if (length(unique(pdat[["Serie"]])) >= 2) {
  pwc.long <- emmeans_test(dplyr::group_by_at(wdat.long, c("grupo","Serie")),
                           score.tde ~ time,
                           p.adjust.method = "bonferroni")
  lpwc[["grupo:Serie"]] <- merge(plyr::rbind.fill(pwc, pwc.long),
                                         lpwc[["grupo:Serie"]],
                                         by=c("grupo","Serie","term",".y.","group1","group2"),
                                         suffixes = c("","'"))
}
grupo Serie term .y. group1 group2 df statistic p p.adj p.adj.signif
Controle 6 ano time score.tde pre pos 2096 1.733 0.083 0.083 ns
Controle 7 ano time score.tde pre pos 2096 2.740 0.006 0.006 **
Controle 8 ano time score.tde pre pos 2096 0.783 0.433 0.433 ns
Controle 9 ano time score.tde pre pos 2096 -0.131 0.896 0.896 ns
Experimental 6 ano time score.tde pre pos 2096 1.773 0.076 0.076 ns
Experimental 7 ano time score.tde pre pos 2096 0.944 0.345 0.345 ns
Experimental 8 ano time score.tde pre pos 2096 -0.962 0.336 0.336 ns
Experimental 9 ano time score.tde pre pos 2096 -1.132 0.258 0.258 ns 
if (length(unique(pdat[["Serie"]])) >= 2) {
  ds <- get.descriptives(wdat, "score.tde.pos", c("grupo","Serie"), covar = "score.tde.pre")
  ds <- merge(ds[ds$variable != "score.tde.pre",],
              ds[ds$variable == "score.tde.pre", !colnames(ds) %in% c("variable")],
              by = c("grupo","Serie"), all.x = T, suffixes = c("", ".score.tde.pre"))
  ds <- merge(get_emmeans(pwcs[["grupo"]]), ds,
              by = c("grupo","Serie"), suffixes = c(".emms", ""))
  ds <- ds[,c("grupo","Serie","n","mean.score.tde.pre","se.score.tde.pre",
              "mean","se","emmean","se.emms","conf.low","conf.high")]
  
  colnames(ds) <- c("grupo","Serie", "N", paste0(c("M","SE")," (pre)"),
                    paste0(c("M","SE"), " (unadj)"),
                    paste0(c("M", "SE"), " (adj)"), "conf.low", "conf.high")
  
  lemms[["grupo:Serie"]] <- merge(ds, lemms[["grupo:Serie"]],
                                          by=c("grupo","Serie"), suffixes = c("","'"))
}
grupo Serie N M (pre) SE (pre) M (unadj) SE (unadj) M (adj) SE (adj) conf.low conf.high
Controle 6 ano 126 28.444 1.766 24.619 1.828 32.399 0.756 30.915 33.882
Controle 7 ano 127 34.024 1.632 28.000 1.773 30.663 0.744 29.204 32.123
Controle 8 ano 86 42.430 1.785 40.337 1.933 35.292 0.906 33.514 37.069
Controle 9 ano 116 46.612 1.611 46.914 1.537 38.034 0.790 36.483 39.584
Experimental 6 ano 150 28.607 1.324 25.020 1.367 32.651 0.694 31.288 34.013
Experimental 7 ano 170 36.294 1.168 34.500 1.284 35.081 0.642 33.822 36.340
Experimental 8 ano 141 39.248 1.525 41.255 1.600 39.128 0.705 37.744 40.512
Experimental 9 ano 140 43.143 1.354 45.514 1.288 39.815 0.713 38.416 41.215

Plots for ancova

if (length(unique(pdat[["Serie"]])) >= 2) {
  ggPlotAoC2(pwcs, "grupo", "Serie", aov, ylab = "Writing (TDE)",
             subtitle = which(aov$Effect == "grupo:Serie"), addParam = "errorbar") +
    ggplot2::scale_color_manual(values = color[["Serie"]])  +
    ggplot2::ylab("Writing (TDE)") +
    theme(axis.title = element_text(size = 14),
          legend.text = element_text(size = 16),
          plot.subtitle = element_text(size = 18)) +
    if (ymin.ci < ymax.ci) ggplot2::ylim(ymin.ci, ymax.ci)
}
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.

if (length(unique(pdat[["Serie"]])) >= 2) {
  ggPlotAoC2(pwcs, "Serie", "grupo", aov, ylab = "Writing (TDE)",
               subtitle = which(aov$Effect == "grupo:Serie"), addParam = "errorbar") +
      ggplot2::scale_color_manual(values = color[["grupo"]])  +
      ggplot2::ylab("Writing (TDE)") +
      theme(axis.title = element_text(size = 14),
            legend.text = element_text(size = 16),
            plot.subtitle = element_text(size = 18)) +
      if (ymin.ci < ymax.ci) ggplot2::ylim(ymin.ci, ymax.ci)
}
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.

if (length(unique(pdat[["Serie"]])) >= 2) {
  plots <- twoWayAncovaBoxPlots(
    wdat, "score.tde.pos", c("grupo","Serie"), aov, pwcs, covar = "score.tde.pre",
    theme = "classic", color = color[["grupo:Serie"]],
    subtitle = which(aov$Effect == "grupo:Serie"))
}
if (length(unique(pdat[["Serie"]])) >= 2) {
  plots[["grupo:Serie"]] + ggplot2::ylab("Writing (TDE)") +
  ggplot2::scale_x_discrete(labels=c('pre', 'pos')) +
  if (ymin < ymax) ggplot2::ylim(ymin, ymax)
}

if (length(unique(pdat[["Serie"]])) >= 2) {
  plots <- twoWayAncovaBoxPlots(
    wdat.long, "score.tde", c("grupo","Serie"), aov, pwc.long,
    pre.post = "time",
    theme = "classic", color = color$prepost)
}
if (length(unique(pdat[["Serie"]])) >= 2) 
  plots[["grupo:Serie"]] + ggplot2::ylab("Writing (TDE)") +
    if (ymin < ymax) ggplot2::ylim(ymin, ymax)

Checking linearity assumption

if (length(unique(pdat[["Serie"]])) >= 2) {
  ggscatter(wdat, x = "score.tde.pre", y = "score.tde.pos", size = 0.5,
            facet.by = c("grupo","Serie"), add = "reg.line")+
    stat_regline_equation(
      aes(label =  paste(..eq.label.., ..rr.label.., sep = "~~~~"))
    ) +
    ggplot2::xlab("Writing (TDE)") +
    ggplot2::ylab("Writing (TDE)") +
    theme(axis.title = element_text(size = 14),
          legend.text = element_text(size = 16),
          plot.subtitle = element_text(size = 18)) +
    if (ymin < ymax) ggplot2::ylim(ymin, ymax)
}

if (length(unique(pdat[["Serie"]])) >= 2) {
  ggscatter(wdat, x = "score.tde.pre", y = "score.tde.pos", size = 0.5,
            color = "grupo", facet.by = "Serie", add = "reg.line")+
    stat_regline_equation(
      aes(label =  paste(..eq.label.., ..rr.label.., sep = "~~~~"), color = grupo)
    ) +
    ggplot2::labs(subtitle = rstatix::get_test_label(aov, detailed = T, row = which(aov$Effect == "grupo:Serie"))) +
    ggplot2::scale_color_manual(values = color[["grupo"]]) +
    ggplot2::xlab("Writing (TDE)") +
    ggplot2::ylab("Writing (TDE)") +
    theme(axis.title = element_text(size = 14),
          legend.text = element_text(size = 16),
          plot.subtitle = element_text(size = 18)) +
    if (ymin < ymax) ggplot2::ylim(ymin, ymax)
}

if (length(unique(pdat[["Serie"]])) >= 2) {
  ggscatter(wdat, x = "score.tde.pre", y = "score.tde.pos", size = 0.5,
            color = "Serie", facet.by = "grupo", add = "reg.line")+
    stat_regline_equation(
      aes(label =  paste(..eq.label.., ..rr.label.., sep = "~~~~"), color = Serie)
    ) +
    ggplot2::labs(subtitle = rstatix::get_test_label(aov, detailed = T, row = which(aov$Effect == "grupo:Serie"))) +
    ggplot2::scale_color_manual(values = color[["Serie"]]) +
    ggplot2::xlab("Writing (TDE)") +
    ggplot2::ylab("Writing (TDE)") +
    theme(axis.title = element_text(size = 14),
          legend.text = element_text(size = 16),
          plot.subtitle = element_text(size = 18)) +
    if (ymin < ymax) ggplot2::ylim(ymin, ymax)
}

Checking normality and homogeneity

if (length(unique(pdat[["Serie"]])) >= 2) 
  res <- augment(lm(score.tde.pos ~ score.tde.pre + grupo*Serie, data = wdat))
if (length(unique(pdat[["Serie"]])) >= 2)
  shapiro_test(res$.resid)
## # A tibble: 1 × 3
##   variable   statistic p.value
##   <chr>          <dbl>   <dbl>
## 1 res$.resid     0.996 0.00727
if (length(unique(pdat[["Serie"]])) >= 2) 
  levene_test(res, .resid ~ grupo*Serie)
## # A tibble: 1 × 4
##     df1   df2 statistic             p
##   <int> <int>     <dbl>         <dbl>
## 1     7  1048      7.87 0.00000000252

ANCOVA and Pairwise for two factors grupo:score.tde.quintile

Without remove non-normal data

pdat = remove_group_data(dat[!is.na(dat[["grupo"]]) & !is.na(dat[["score.tde.quintile"]]),],
                         "score.tde.pos", c("grupo","score.tde.quintile"))
pdat = pdat[pdat[["score.tde.quintile"]] %in% do.call(
  intersect, lapply(unique(pdat[["grupo"]]), FUN = function(x) {
    unique(pdat[["score.tde.quintile"]][which(pdat[["grupo"]] == x)])
  })),]
pdat[["grupo"]] = factor(pdat[["grupo"]], level[["grupo"]])
pdat[["score.tde.quintile"]] = factor(
  pdat[["score.tde.quintile"]],
  level[["score.tde.quintile"]][level[["score.tde.quintile"]] %in% unique(pdat[["score.tde.quintile"]])])

pdat.long <- rbind(pdat[,c("id","grupo","score.tde.quintile")], pdat[,c("id","grupo","score.tde.quintile")])
pdat.long[["time"]] <- c(rep("pre", nrow(pdat)), rep("pos", nrow(pdat)))
pdat.long[["time"]] <- factor(pdat.long[["time"]], c("pre","pos"))
pdat.long[["score.tde"]] <- c(pdat[["score.tde.pre"]], pdat[["score.tde.pos"]])

if (length(unique(pdat[["score.tde.quintile"]])) >= 2) {
  aov = anova_test(pdat, score.tde.pos ~ score.tde.pre + grupo*score.tde.quintile)
  laov[["grupo:score.tde.quintile"]] <- get_anova_table(aov)
}
if (length(unique(pdat[["score.tde.quintile"]])) >= 2) {
  pwcs <- list()
  pwcs[["score.tde.quintile"]] <- emmeans_test(
    group_by(pdat, grupo), score.tde.pos ~ score.tde.quintile,
    covariate = score.tde.pre, p.adjust.method = "bonferroni")
  pwcs[["grupo"]] <- emmeans_test(
    group_by(pdat, score.tde.quintile), score.tde.pos ~ grupo,
    covariate = score.tde.pre, p.adjust.method = "bonferroni")
  
  pwc <- plyr::rbind.fill(pwcs[["grupo"]], pwcs[["score.tde.quintile"]])
  pwc <- pwc[,c("grupo","score.tde.quintile", colnames(pwc)[!colnames(pwc) %in% c("grupo","score.tde.quintile")])]
}
if (length(unique(pdat[["score.tde.quintile"]])) >= 2) {
  pwc.long <- emmeans_test(dplyr::group_by_at(pdat.long, c("grupo","score.tde.quintile")),
                           score.tde ~ time,
                           p.adjust.method = "bonferroni")
  lpwc[["grupo:score.tde.quintile"]] <- plyr::rbind.fill(pwc, pwc.long)
}
if (length(unique(pdat[["score.tde.quintile"]])) >= 2) {
  ds <- get.descriptives(pdat, "score.tde.pos", c("grupo","score.tde.quintile"), covar = "score.tde.pre")
  ds <- merge(ds[ds$variable != "score.tde.pre",],
              ds[ds$variable == "score.tde.pre", !colnames(ds) %in% c("variable")],
              by = c("grupo","score.tde.quintile"), all.x = T, suffixes = c("", ".score.tde.pre"))
  ds <- merge(get_emmeans(pwcs[["grupo"]]), ds,
              by = c("grupo","score.tde.quintile"), suffixes = c(".emms", ""))
  ds <- ds[,c("grupo","score.tde.quintile","n","mean.score.tde.pre","se.score.tde.pre","mean","se",
              "emmean","se.emms","conf.low","conf.high")]
  
  colnames(ds) <- c("grupo","score.tde.quintile", "N", paste0(c("M","SE")," (pre)"),
                    paste0(c("M","SE"), " (unadj)"),
                    paste0(c("M", "SE"), " (adj)"), "conf.low", "conf.high")
  
  lemms[["grupo:score.tde.quintile"]] <- ds
}

Computing ANCOVA and PairWise After removing non-normal data (OK)

if (length(unique(pdat[["score.tde.quintile"]])) >= 2) {
  wdat = pdat 
  
  res = residuals(lm(score.tde.pos ~ score.tde.pre + grupo*score.tde.quintile, data = wdat))
  non.normal = getNonNormal(res, wdat$id, plimit = 0.05)
  
  wdat = wdat[!wdat$id %in% non.normal,]
  
  wdat.long <- rbind(wdat[,c("id","grupo","score.tde.quintile")], wdat[,c("id","grupo","score.tde.quintile")])
  wdat.long[["time"]] <- c(rep("pre", nrow(wdat)), rep("pos", nrow(wdat)))
  wdat.long[["time"]] <- factor(wdat.long[["time"]], c("pre","pos"))
  wdat.long[["score.tde"]] <- c(wdat[["score.tde.pre"]], wdat[["score.tde.pos"]])
  
  
  ldat[["grupo:score.tde.quintile"]] = wdat
  
  (non.normal)
}
##   [1] "P984"  "P1128" "P962"  "P1129" "P3637" "P1971" "P1117" "P1126" "P2950"
##  [10] "P2883" "P2910" "P3019" "P2982" "P2904" "P2864" "P2959" "P3054" "P1139"
##  [19] "P2835" "P921"  "P2983" "P2994" "P2880" "P2917" "P2848" "P2866" "P3005"
##  [28] "P929"  "P2870" "P572"  "P3674" "P3548" "P3021" "P809"  "P3080" "P2997"
##  [37] "P3666" "P3660" "P2867" "P1878" "P2858" "P2846" "P1840" "P3007" "P3545"
##  [46] "P3533" "P2861" "P2993" "P2953" "P1111" "P2969" "P1804" "P2876" "P1018"
##  [55] "P3228" "P2854" "P2871" "P2929" "P908"  "P2964" "P1046" "P2190" "P2879"
##  [64] "P2948" "P1118" "P2967" "P2956" "P2973" "P2946" "P2865" "P2831" "P3000"
##  [73] "P2995" "P976"  "P2937" "P1145" "P3029" "P2971" "P2951" "P2886" "P2913"
##  [82] "P2977" "P1885" "P3026" "P931"  "P3475" "P924"  "P2952" "P2974" "P1024"
##  [91] "P3520" "P3136" "P914"  "P3280" "P1130" "P1825" "P2850" "P2888" "P1968"
## [100] "P1083" "P1887" "P2905" "P3459" "P3574"
if (length(unique(pdat[["score.tde.quintile"]])) >= 2) {
  aov = anova_test(wdat, score.tde.pos ~ score.tde.pre + grupo*score.tde.quintile)
  laov[["grupo:score.tde.quintile"]] <- merge(get_anova_table(aov), laov[["grupo:score.tde.quintile"]],
                                         by="Effect", suffixes = c("","'"))
  df = get_anova_table(aov)
}
Effect DFn DFd F p p<.05 ges
score.tde.pre 1 1006 264.709 0.000 * 0.208
grupo 1 1006 37.367 0.000 * 0.036
score.tde.quintile 4 1006 3.315 0.010 * 0.013
grupo:score.tde.quintile 4 1006 2.738 0.028 * 0.011 
if (length(unique(pdat[["score.tde.quintile"]])) >= 2) {
  pwcs <- list()
  pwcs[["score.tde.quintile"]] <- emmeans_test(
    group_by(wdat, grupo), score.tde.pos ~ score.tde.quintile,
    covariate = score.tde.pre, p.adjust.method = "bonferroni")
  pwcs[["grupo"]] <- emmeans_test(
    group_by(wdat, score.tde.quintile), score.tde.pos ~ grupo,
    covariate = score.tde.pre, p.adjust.method = "bonferroni")
  
  pwc <- plyr::rbind.fill(pwcs[["grupo"]], pwcs[["score.tde.quintile"]])
  pwc <- pwc[,c("grupo","score.tde.quintile", colnames(pwc)[!colnames(pwc) %in% c("grupo","score.tde.quintile")])]
}
grupo score.tde.quintile term .y. group1 group2 df statistic p p.adj p.adj.signif
1st quintile score.tde.pre*grupo score.tde.pos Controle Experimental 1006 -2.176 0.030 0.030 *
2nd quintile score.tde.pre*grupo score.tde.pos Controle Experimental 1006 -2.614 0.009 0.009 **
3rd quintile score.tde.pre*grupo score.tde.pos Controle Experimental 1006 -4.957 0.000 0.000 ****
4th quintile score.tde.pre*grupo score.tde.pos Controle Experimental 1006 -3.161 0.002 0.002 **
5th quintile score.tde.pre*grupo score.tde.pos Controle Experimental 1006 -1.453 0.147 0.147 ns
Controle score.tde.pre*score.tde.quintile score.tde.pos 1st quintile 2nd quintile 1006 1.565 0.118 1.000 ns
Controle score.tde.pre*score.tde.quintile score.tde.pos 1st quintile 3rd quintile 1006 2.102 0.036 0.358 ns
Controle score.tde.pre*score.tde.quintile score.tde.pos 1st quintile 4th quintile 1006 -0.118 0.906 1.000 ns
Controle score.tde.pre*score.tde.quintile score.tde.pos 1st quintile 5th quintile 1006 -0.370 0.711 1.000 ns
Controle score.tde.pre*score.tde.quintile score.tde.pos 2nd quintile 3rd quintile 1006 1.201 0.230 1.000 ns
Controle score.tde.pre*score.tde.quintile score.tde.pos 2nd quintile 4th quintile 1006 -1.486 0.138 1.000 ns
Controle score.tde.pre*score.tde.quintile score.tde.pos 2nd quintile 5th quintile 1006 -1.515 0.130 1.000 ns
Controle score.tde.pre*score.tde.quintile score.tde.pos 3rd quintile 4th quintile 1006 -3.359 0.001 0.008 **
Controle score.tde.pre*score.tde.quintile score.tde.pos 3rd quintile 5th quintile 1006 -3.197 0.001 0.014 *
Controle score.tde.pre*score.tde.quintile score.tde.pos 4th quintile 5th quintile 1006 -0.703 0.482 1.000 ns
Experimental score.tde.pre*score.tde.quintile score.tde.pos 1st quintile 2nd quintile 1006 0.882 0.378 1.000 ns
Experimental score.tde.pre*score.tde.quintile score.tde.pos 1st quintile 3rd quintile 1006 -0.005 0.996 1.000 ns
Experimental score.tde.pre*score.tde.quintile score.tde.pos 1st quintile 4th quintile 1006 -0.482 0.630 1.000 ns
Experimental score.tde.pre*score.tde.quintile score.tde.pos 1st quintile 5th quintile 1006 -0.099 0.921 1.000 ns
Experimental score.tde.pre*score.tde.quintile score.tde.pos 2nd quintile 3rd quintile 1006 -1.014 0.311 1.000 ns
Experimental score.tde.pre*score.tde.quintile score.tde.pos 2nd quintile 4th quintile 1006 -1.503 0.133 1.000 ns
Experimental score.tde.pre*score.tde.quintile score.tde.pos 2nd quintile 5th quintile 1006 -0.722 0.470 1.000 ns
Experimental score.tde.pre*score.tde.quintile score.tde.pos 3rd quintile 4th quintile 1006 -0.991 0.322 1.000 ns
Experimental score.tde.pre*score.tde.quintile score.tde.pos 3rd quintile 5th quintile 1006 -0.183 0.855 1.000 ns
Experimental score.tde.pre*score.tde.quintile score.tde.pos 4th quintile 5th quintile 1006 0.707 0.480 1.000 ns 
if (length(unique(pdat[["score.tde.quintile"]])) >= 2) {
  pwc.long <- emmeans_test(dplyr::group_by_at(wdat.long, c("grupo","score.tde.quintile")),
                           score.tde ~ time,
                           p.adjust.method = "bonferroni")
  lpwc[["grupo:score.tde.quintile"]] <- merge(plyr::rbind.fill(pwc, pwc.long),
                                         lpwc[["grupo:score.tde.quintile"]],
                                         by=c("grupo","score.tde.quintile","term",".y.","group1","group2"),
                                         suffixes = c("","'"))
}
grupo score.tde.quintile term .y. group1 group2 df statistic p p.adj p.adj.signif
Controle 1st quintile time score.tde pre pos 2014 1.024 0.306 0.306 ns
Controle 2nd quintile time score.tde pre pos 2014 3.105 0.002 0.002 **
Controle 3rd quintile time score.tde pre pos 2014 4.794 0.000 0.000 ****
Controle 4th quintile time score.tde pre pos 2014 3.089 0.002 0.002 **
Controle 5th quintile time score.tde pre pos 2014 2.704 0.007 0.007 **
Experimental 1st quintile time score.tde pre pos 2014 -1.416 0.157 0.157 ns
Experimental 2nd quintile time score.tde pre pos 2014 0.824 0.410 0.410 ns
Experimental 3rd quintile time score.tde pre pos 2014 0.152 0.879 0.879 ns
Experimental 4th quintile time score.tde pre pos 2014 -0.569 0.569 0.569 ns
Experimental 5th quintile time score.tde pre pos 2014 1.120 0.263 0.263 ns 
if (length(unique(pdat[["score.tde.quintile"]])) >= 2) {
  ds <- get.descriptives(wdat, "score.tde.pos", c("grupo","score.tde.quintile"), covar = "score.tde.pre")
  ds <- merge(ds[ds$variable != "score.tde.pre",],
              ds[ds$variable == "score.tde.pre", !colnames(ds) %in% c("variable")],
              by = c("grupo","score.tde.quintile"), all.x = T, suffixes = c("", ".score.tde.pre"))
  ds <- merge(get_emmeans(pwcs[["grupo"]]), ds,
              by = c("grupo","score.tde.quintile"), suffixes = c(".emms", ""))
  ds <- ds[,c("grupo","score.tde.quintile","n","mean.score.tde.pre","se.score.tde.pre",
              "mean","se","emmean","se.emms","conf.low","conf.high")]
  
  colnames(ds) <- c("grupo","score.tde.quintile", "N", paste0(c("M","SE")," (pre)"),
                    paste0(c("M","SE"), " (unadj)"),
                    paste0(c("M", "SE"), " (adj)"), "conf.low", "conf.high")
  
  lemms[["grupo:score.tde.quintile"]] <- merge(ds, lemms[["grupo:score.tde.quintile"]],
                                          by=c("grupo","score.tde.quintile"), suffixes = c("","'"))
}
grupo score.tde.quintile N M (pre) SE (pre) M (unadj) SE (unadj) M (adj) SE (adj) conf.low conf.high
Controle 1st quintile 112 9.152 0.522 8.205 0.750 34.585 1.783 31.086 38.083
Controle 2nd quintile 51 24.392 0.530 20.137 1.288 32.089 1.322 29.495 34.682
Controle 3rd quintile 42 38.905 0.489 31.667 1.800 29.879 1.216 27.493 32.265
Controle 4th quintile 125 47.544 0.223 44.840 0.850 34.874 0.932 33.046 36.702
Controle 5th quintile 116 58.819 0.505 56.362 0.792 35.722 1.463 32.851 38.593
Experimental 1st quintile 108 9.602 0.522 10.935 0.788 36.889 1.765 33.425 40.352
Experimental 2nd quintile 101 25.842 0.347 25.040 0.976 35.619 1.016 33.625 37.613
Experimental 3rd quintile 116 37.371 0.301 37.233 0.859 36.898 0.729 35.467 38.328
Experimental 4th quintile 120 46.800 0.230 47.308 0.725 38.047 0.915 36.251 39.842
Experimental 5th quintile 126 58.563 0.448 57.587 0.700 37.189 1.435 34.372 40.006

Plots for ancova

if (length(unique(pdat[["score.tde.quintile"]])) >= 2) {
  ggPlotAoC2(pwcs, "grupo", "score.tde.quintile", aov, ylab = "Writing (TDE)",
             subtitle = which(aov$Effect == "grupo:score.tde.quintile"), addParam = "errorbar") +
    ggplot2::scale_color_manual(values = color[["score.tde.quintile"]])  +
    ggplot2::ylab("Writing (TDE)") +
    theme(axis.title = element_text(size = 14),
          legend.text = element_text(size = 16),
          plot.subtitle = element_text(size = 18)) +
    if (ymin.ci < ymax.ci) ggplot2::ylim(ymin.ci, ymax.ci)
}
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.

if (length(unique(pdat[["score.tde.quintile"]])) >= 2) {
  ggPlotAoC2(pwcs, "score.tde.quintile", "grupo", aov, ylab = "Writing (TDE)",
               subtitle = which(aov$Effect == "grupo:score.tde.quintile"), addParam = "errorbar") +
      ggplot2::scale_color_manual(values = color[["grupo"]])  +
      ggplot2::ylab("Writing (TDE)") +
      theme(axis.title = element_text(size = 14),
            legend.text = element_text(size = 16),
            plot.subtitle = element_text(size = 18)) +
      if (ymin.ci < ymax.ci) ggplot2::ylim(ymin.ci, ymax.ci)
}
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.

if (length(unique(pdat[["score.tde.quintile"]])) >= 2) {
  plots <- twoWayAncovaBoxPlots(
    wdat, "score.tde.pos", c("grupo","score.tde.quintile"), aov, pwcs, covar = "score.tde.pre",
    theme = "classic", color = color[["grupo:score.tde.quintile"]],
    subtitle = which(aov$Effect == "grupo:score.tde.quintile"))
}
if (length(unique(pdat[["score.tde.quintile"]])) >= 2) {
  plots[["grupo:score.tde.quintile"]] + ggplot2::ylab("Writing (TDE)") +
  ggplot2::scale_x_discrete(labels=c('pre', 'pos')) +
  if (ymin < ymax) ggplot2::ylim(ymin, ymax)
}
## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's colour values.

if (length(unique(pdat[["score.tde.quintile"]])) >= 2) {
  plots <- twoWayAncovaBoxPlots(
    wdat.long, "score.tde", c("grupo","score.tde.quintile"), aov, pwc.long,
    pre.post = "time",
    theme = "classic", color = color$prepost)
}
if (length(unique(pdat[["score.tde.quintile"]])) >= 2) 
  plots[["grupo:score.tde.quintile"]] + ggplot2::ylab("Writing (TDE)") +
    if (ymin < ymax) ggplot2::ylim(ymin, ymax)

Checking linearity assumption

if (length(unique(pdat[["score.tde.quintile"]])) >= 2) {
  ggscatter(wdat, x = "score.tde.pre", y = "score.tde.pos", size = 0.5,
            facet.by = c("grupo","score.tde.quintile"), add = "reg.line")+
    stat_regline_equation(
      aes(label =  paste(..eq.label.., ..rr.label.., sep = "~~~~"))
    ) +
    ggplot2::xlab("Writing (TDE)") +
    ggplot2::ylab("Writing (TDE)") +
    theme(axis.title = element_text(size = 14),
          legend.text = element_text(size = 16),
          plot.subtitle = element_text(size = 18)) +
    if (ymin < ymax) ggplot2::ylim(ymin, ymax)
}

if (length(unique(pdat[["score.tde.quintile"]])) >= 2) {
  ggscatter(wdat, x = "score.tde.pre", y = "score.tde.pos", size = 0.5,
            color = "grupo", facet.by = "score.tde.quintile", add = "reg.line")+
    stat_regline_equation(
      aes(label =  paste(..eq.label.., ..rr.label.., sep = "~~~~"), color = grupo)
    ) +
    ggplot2::labs(subtitle = rstatix::get_test_label(aov, detailed = T, row = which(aov$Effect == "grupo:score.tde.quintile"))) +
    ggplot2::scale_color_manual(values = color[["grupo"]]) +
    ggplot2::xlab("Writing (TDE)") +
    ggplot2::ylab("Writing (TDE)") +
    theme(axis.title = element_text(size = 14),
          legend.text = element_text(size = 16),
          plot.subtitle = element_text(size = 18)) +
    if (ymin < ymax) ggplot2::ylim(ymin, ymax)
}

if (length(unique(pdat[["score.tde.quintile"]])) >= 2) {
  ggscatter(wdat, x = "score.tde.pre", y = "score.tde.pos", size = 0.5,
            color = "score.tde.quintile", facet.by = "grupo", add = "reg.line")+
    stat_regline_equation(
      aes(label =  paste(..eq.label.., ..rr.label.., sep = "~~~~"), color = score.tde.quintile)
    ) +
    ggplot2::labs(subtitle = rstatix::get_test_label(aov, detailed = T, row = which(aov$Effect == "grupo:score.tde.quintile"))) +
    ggplot2::scale_color_manual(values = color[["score.tde.quintile"]]) +
    ggplot2::xlab("Writing (TDE)") +
    ggplot2::ylab("Writing (TDE)") +
    theme(axis.title = element_text(size = 14),
          legend.text = element_text(size = 16),
          plot.subtitle = element_text(size = 18)) +
    if (ymin < ymax) ggplot2::ylim(ymin, ymax)
}

Checking normality and homogeneity

if (length(unique(pdat[["score.tde.quintile"]])) >= 2) 
  res <- augment(lm(score.tde.pos ~ score.tde.pre + grupo*score.tde.quintile, data = wdat))
if (length(unique(pdat[["score.tde.quintile"]])) >= 2)
  shapiro_test(res$.resid)
## # A tibble: 1 × 3
##   variable   statistic p.value
##   <chr>          <dbl>   <dbl>
## 1 res$.resid     0.996  0.0129
if (length(unique(pdat[["score.tde.quintile"]])) >= 2) 
  levene_test(res, .resid ~ grupo*score.tde.quintile)
## # A tibble: 1 × 4
##     df1   df2 statistic           p
##   <int> <int>     <dbl>       <dbl>
## 1     9  1007      5.39 0.000000314

Summary of Results

Descriptive Statistics

df <- get.descriptives(ldat[["grupo"]], c(dv.pre, dv.pos), c("grupo"), 
                       include.global = T, symmetry.test = T, normality.test = F)
df <- plyr::rbind.fill(
  df, do.call(plyr::rbind.fill, lapply(lfatores2, FUN = function(f) {
    if (nrow(dat) > 0 && sum(!is.na(unique(dat[[f]]))) > 1 && paste0("grupo:",f) %in% names(ldat))
      get.descriptives(ldat[[paste0("grupo:",f)]], c(dv.pre,dv.pos), c("grupo", f),
                       symmetry.test = T, normality.test = F)
    }))
)
df <- df[,c(fatores1[fatores1 %in% colnames(df)],"variable",
             colnames(df)[!colnames(df) %in% c(fatores1,"variable")])]
grupo Sexo Zona Cor.Raca Serie score.tde.quintile variable n mean median min max sd se ci iqr symmetry skewness kurtosis
Controle score.tde.pre 424 37.351 45.0 0 75 20.069 0.975 1.916 34.00 YES -0.404 -1.165
Experimental score.tde.pre 557 36.770 39.0 0 73 17.478 0.741 1.455 26.00 YES -0.313 -0.796
score.tde.pre 981 37.021 42.0 0 75 18.635 0.595 1.168 29.00 YES -0.358 -0.971
Controle score.tde.pos 424 35.432 40.0 0 73 20.965 1.018 2.001 38.00 YES -0.277 -1.217
Experimental score.tde.pos 557 37.230 40.0 0 74 18.335 0.777 1.526 27.00 YES -0.294 -0.858
score.tde.pos 981 36.453 40.0 0 74 19.525 0.623 1.223 31.00 YES -0.305 -1.008
Controle F score.tde.pre 208 40.635 47.0 0 72 18.428 1.278 2.519 28.00 NO -0.660 -0.764
Controle M score.tde.pre 211 33.934 42.0 0 75 21.161 1.457 2.872 38.00 YES -0.126 -1.375
Experimental F score.tde.pre 276 39.623 43.5 0 73 17.500 1.053 2.074 26.00 YES -0.490 -0.692
Experimental M score.tde.pre 279 34.082 36.0 0 71 17.200 1.030 2.027 26.00 YES -0.174 -0.823
Controle F score.tde.pos 208 38.615 44.0 0 71 19.482 1.351 2.663 30.00 NO -0.532 -0.975
Controle M score.tde.pos 211 32.393 33.0 0 73 22.040 1.517 2.991 38.00 YES -0.018 -1.332
Experimental F score.tde.pos 276 40.272 45.0 0 73 18.658 1.123 2.211 28.25 YES -0.484 -0.805
Experimental M score.tde.pos 279 34.484 37.0 0 74 17.582 1.053 2.072 25.00 YES -0.205 -0.757
Controle Rural score.tde.pre 228 35.956 43.5 0 69 20.084 1.330 2.621 37.00 YES -0.398 -1.296
Controle Urbana score.tde.pre 101 39.743 45.0 1 75 18.324 1.823 3.617 28.00 YES -0.349 -0.771
Experimental Rural score.tde.pre 260 34.215 34.5 0 73 17.945 1.113 2.192 27.50 YES -0.096 -0.905
Experimental Urbana score.tde.pre 151 38.318 42.0 0 71 17.068 1.389 2.745 23.50 NO -0.526 -0.563
Controle Rural score.tde.pos 228 34.978 42.0 0 72 21.528 1.426 2.809 38.25 YES -0.286 -1.280
Controle Urbana score.tde.pos 101 31.099 30.0 0 73 20.538 2.044 4.054 36.00 YES 0.126 -1.232
Experimental Rural score.tde.pos 260 34.777 36.0 0 72 19.125 1.186 2.336 32.00 YES -0.055 -1.056
Experimental Urbana score.tde.pos 151 36.099 40.0 0 74 18.449 1.501 2.967 25.50 YES -0.391 -0.821
Controle Parda score.tde.pre 150 36.320 43.0 0 66 19.333 1.578 3.119 33.00 NO -0.563 -1.085
Controle Indígena score.tde.pre 11 42.000 46.0 4 65 17.297 5.215 11.621 9.50 NO -1.005 -0.104
Controle Branca score.tde.pre 45 41.378 46.0 2 67 17.619 2.626 5.293 15.00 NO -0.795 -0.473
Experimental Parda score.tde.pre 175 33.851 35.0 0 68 17.760 1.343 2.650 28.00 YES -0.110 -0.955
Experimental Indígena score.tde.pre 15 27.067 26.0 6 54 17.065 4.406 9.450 27.00 YES 0.126 -1.422
Experimental Branca score.tde.pre 58 34.517 35.5 0 73 20.186 2.651 5.308 34.50 YES -0.030 -1.118
Controle Parda score.tde.pos 150 33.887 39.5 0 68 20.756 1.695 3.349 39.50 YES -0.286 -1.375
Controle Indígena score.tde.pos 11 45.636 52.0 4 63 19.931 6.009 13.390 19.50 NO -0.968 -0.645
Controle Branca score.tde.pos 45 36.600 43.0 0 72 20.175 3.007 6.061 28.00 YES -0.398 -0.997
Experimental Parda score.tde.pos 175 33.514 32.0 0 71 19.317 1.460 2.882 31.00 YES 0.015 -1.092
Experimental Indígena score.tde.pos 15 30.333 26.0 1 57 17.715 4.574 9.810 30.50 YES 0.117 -1.542
Experimental Branca score.tde.pos 58 34.586 34.5 0 72 20.910 2.746 5.498 32.00 YES 0.002 -1.163
Controle 6 ano score.tde.pre 126 28.444 25.0 0 63 19.827 1.766 3.496 37.00 YES 0.114 -1.549
Controle 7 ano score.tde.pre 127 34.024 41.0 0 69 18.393 1.632 3.230 30.00 YES -0.192 -1.277
Controle 8 ano score.tde.pre 86 42.430 46.5 0 72 16.556 1.785 3.550 18.50 NO -0.821 0.143
Controle 9 ano score.tde.pre 116 46.612 51.0 0 75 17.352 1.611 3.191 14.25 NO -1.081 0.454
Experimental 6 ano score.tde.pre 150 28.607 29.0 0 65 16.217 1.324 2.617 24.75 YES 0.006 -0.955
Experimental 7 ano score.tde.pre 170 36.294 37.0 0 68 15.233 1.168 2.306 22.00 YES -0.132 -0.765
Experimental 8 ano score.tde.pre 141 39.248 44.0 0 73 18.114 1.525 3.016 26.00 NO -0.545 -0.691
Experimental 9 ano score.tde.pre 140 43.143 45.0 0 71 16.016 1.354 2.676 22.25 NO -0.646 -0.066
Controle 6 ano score.tde.pos 126 24.619 20.0 0 66 20.516 1.828 3.617 40.75 YES 0.335 -1.356
Controle 7 ano score.tde.pos 127 28.000 24.0 0 70 19.979 1.773 3.508 34.50 YES 0.321 -1.224
Controle 8 ano score.tde.pos 86 40.337 43.5 0 71 17.930 1.933 3.844 26.75 NO -0.595 -0.451
Controle 9 ano score.tde.pos 116 46.914 50.0 0 73 16.558 1.537 3.045 18.25 NO -0.869 0.170
Experimental 6 ano score.tde.pos 150 25.020 23.0 0 65 16.738 1.367 2.701 27.00 YES 0.452 -0.667
Experimental 7 ano score.tde.pos 170 34.500 37.0 0 67 16.741 1.284 2.535 25.75 YES -0.209 -0.863
Experimental 8 ano score.tde.pos 141 41.255 45.0 0 73 18.995 1.600 3.163 26.00 NO -0.561 -0.642
Experimental 9 ano score.tde.pos 140 45.514 49.0 0 74 15.245 1.288 2.547 19.50 NO -0.761 0.330
Controle 1st quintile score.tde.pre 112 9.152 10.0 0 18 5.522 0.522 1.034 8.50 YES -0.089 -1.225
Controle 2nd quintile score.tde.pre 51 24.392 25.0 19 31 3.785 0.530 1.064 6.50 YES -0.049 -1.183
Controle 3rd quintile score.tde.pre 42 38.905 40.0 32 42 3.169 0.489 0.987 5.75 NO -0.602 -1.073
Controle 4th quintile score.tde.pre 125 47.544 48.0 43 51 2.497 0.223 0.442 4.00 YES -0.303 -0.977
Controle 5th quintile score.tde.pre 116 58.819 58.0 52 75 5.437 0.505 1.000 8.00 NO 0.749 -0.080
Experimental 1st quintile score.tde.pre 108 9.602 10.0 0 18 5.426 0.522 1.035 10.00 YES -0.174 -1.164
Experimental 2nd quintile score.tde.pre 101 25.842 26.0 19 31 3.489 0.347 0.689 6.00 YES -0.241 -1.059
Experimental 3rd quintile score.tde.pre 116 37.371 38.0 32 42 3.245 0.301 0.597 6.00 YES -0.131 -1.302
Experimental 4th quintile score.tde.pre 120 46.800 47.0 43 51 2.523 0.230 0.456 4.00 YES 0.068 -1.184
Experimental 5th quintile score.tde.pre 126 58.563 58.0 52 73 5.033 0.448 0.887 8.00 NO 0.662 -0.409
Controle 1st quintile score.tde.pos 112 8.205 5.0 0 28 7.936 0.750 1.486 12.25 NO 0.705 -0.691
Controle 2nd quintile score.tde.pos 51 20.137 20.0 4 45 9.200 1.288 2.588 11.50 YES 0.436 0.070
Controle 3rd quintile score.tde.pos 42 31.667 30.0 11 59 11.665 1.800 3.635 17.75 YES 0.319 -0.776
Controle 4th quintile score.tde.pos 125 44.840 46.0 18 63 9.506 0.850 1.683 14.00 NO -0.520 -0.442
Controle 5th quintile score.tde.pos 116 56.362 56.5 31 73 8.531 0.792 1.569 10.25 YES -0.360 -0.038
Experimental 1st quintile score.tde.pos 108 10.935 10.0 0 36 8.185 0.788 1.561 12.00 NO 0.643 -0.292
Experimental 2nd quintile score.tde.pos 101 25.040 25.0 4 47 9.804 0.976 1.935 14.00 YES -0.110 -0.678
Experimental 3rd quintile score.tde.pos 116 37.233 39.0 16 55 9.256 0.859 1.702 14.25 YES -0.215 -0.813
Experimental 4th quintile score.tde.pos 120 47.308 48.0 22 69 7.944 0.725 1.436 9.25 YES -0.124 0.510
Experimental 5th quintile score.tde.pos 126 57.587 57.0 39 74 7.862 0.700 1.386 11.75 YES -0.167 -0.553

ANCOVA Table Comparison

df <- do.call(plyr::rbind.fill, laov)
df <- df[!duplicated(df$Effect),]
Effect DFn DFd F p p<.05 ges DFn’ DFd’ F’ p’ p<.05’ ges’
1 grupo 1 978 24.670 0.000 * 0.025 1 1118 13.989 0.000 * 0.012
2 score.tde.pre 1 978 5874.442 0.000 * 0.857 1 1118 2814.834 0.000 * 0.716
4 grupo:Sexo 1 969 0.548 0.459 0.001 1 1116 0.305 0.581 0.000
6 Sexo 1 969 0.058 0.810 0.000 1 1116 0.148 0.700 0.000
8 grupo:Zona 1 735 14.155 0.000 * 0.019 1 798 9.749 0.002 * 0.012
10 Zona 1 735 54.679 0.000 * 0.069 1 798 35.734 0.000 * 0.043
11 Cor.Raca 2 447 3.357 0.036 * 0.015 2 482 3.984 0.019 * 0.016
13 grupo:Cor.Raca 2 447 1.174 0.310 0.005 2 482 1.222 0.296 0.005
16 grupo:Serie 3 1047 3.528 0.015 * 0.010 3 1112 2.801 0.039 * 0.008
18 Serie 3 1047 34.727 0.000 * 0.090 3 1112 30.848 0.000 * 0.077
20 grupo:score.tde.quintile 4 1006 2.738 0.028 * 0.011 4 1110 1.400 0.232 0.005
22 score.tde.quintile 4 1006 3.315 0.010 * 0.013 4 1110 3.862 0.004 * 0.014

PairWise Table Comparison

df <- do.call(plyr::rbind.fill, lpwc)
df <- df[,c(names(lfatores)[names(lfatores) %in% colnames(df)],
            names(df)[!names(df) %in% c(names(lfatores),"term",".y.")])]
grupo Sexo Zona Cor.Raca Serie score.tde.quintile group1 group2 df statistic p p.adj p.adj.signif df’ statistic’ p’ p.adj’ p.adj.signif’
Controle pre pos 1958 1.465 0.143 0.143 ns 2238 3.143 0.002 0.002 **
Experimental pre pos 1958 -0.402 0.688 0.688 ns 2238 1.295 0.195 0.195 ns
Controle Experimental 978 -4.967 0.000 0.000 **** 1118 -3.740 0.000 0.000 ***
Controle F pre pos 1940 1.089 0.276 0.276 ns 2234 2.427 0.015 0.015 *
Controle M pre pos 1940 0.837 0.403 0.403 ns 2234 2.061 0.039 0.039 *
Controle F M 969 -0.398 0.691 0.691 ns 1116 -0.159 0.874 0.874 ns
Experimental F pre pos 1940 -0.403 0.687 0.687 ns 2234 0.857 0.391 0.391 ns
Experimental M pre pos 1940 -0.251 0.802 0.802 ns 2234 0.994 0.320 0.320 ns
Experimental F M 969 0.663 0.507 0.507 ns 1116 0.652 0.514 0.514 ns
F Controle Experimental 969 -3.977 0.000 0.000 **** 1116 -3.048 0.002 0.002 **
M Controle Experimental 969 -2.953 0.003 0.003 ** 1116 -2.236 0.026 0.026 *
Controle Rural Urbana 735 7.675 0.000 0.000 **** 798 6.249 0.000 0.000 ****
Controle Rural pre pos 1472 0.543 0.588 0.588 ns 1598 1.120 0.263 0.263 ns
Controle Urbana pre pos 1472 3.191 0.001 0.001 ** 1598 3.728 0.000 0.000 ***
Experimental Rural Urbana 735 3.188 0.001 0.001 ** 798 2.559 0.011 0.011 *
Experimental Rural pre pos 1472 -0.333 0.739 0.739 ns 1598 0.603 0.547 0.547 ns
Experimental Urbana pre pos 1472 1.001 0.317 0.317 ns 1598 1.793 0.073 0.073 ns
Rural Controle Experimental 735 -2.009 0.045 0.045 * 798 -0.977 0.329 0.329 ns
Urbana Controle Experimental 735 -6.023 0.000 0.000 **** 798 -4.521 0.000 0.000 ****
Controle Branca pre pos 896 1.172 0.241 0.241 ns 966 1.861 0.063 0.063 ns
Controle Indígena pre pos 896 -0.441 0.659 0.659 ns 966 -0.443 0.658 0.658 ns
Controle Indígena Branca 447 2.687 0.007 0.022 * 482 2.857 0.004 0.013 *
Controle Parda Branca 447 1.334 0.183 0.549 ns 482 1.348 0.178 0.535 ns
Controle Parda Indígena 447 -2.163 0.031 0.093 ns 482 -2.350 0.019 0.057 ns
Controle Parda pre pos 896 1.090 0.276 0.276 ns 966 2.002 0.046 0.046 *
Experimental Branca pre pos 896 -0.019 0.985 0.985 ns 966 0.377 0.707 0.707 ns
Experimental Indígena pre pos 896 -0.463 0.644 0.644 ns 966 -0.465 0.642 0.642 ns
Experimental Indígena Branca 447 1.060 0.290 0.870 ns 482 1.201 0.230 0.691 ns
Experimental Parda Branca 447 -0.308 0.759 1.000 ns 482 -0.313 0.755 1.000 ns
Experimental Parda Indígena 447 -1.314 0.190 0.569 ns 482 -1.461 0.145 0.434 ns
Experimental Parda pre pos 896 0.163 0.870 0.870 ns 966 0.915 0.360 0.360 ns
Branca Controle Experimental 447 -2.442 0.015 0.015 * 482 -2.420 0.016 0.016 *
Indígena Controle Experimental 447 0.274 0.784 0.784 ns 482 0.354 0.723 0.723 ns
Parda Controle Experimental 447 -1.911 0.057 0.057 ns 482 -1.857 0.064 0.064 ns
Controle 6 ano pre pos 2096 1.733 0.083 0.083 ns 2226 2.226 0.026 0.026 *
Controle 7 ano pre pos 2096 2.740 0.006 0.006 ** 2226 3.047 0.002 0.002 **
Controle 8 ano pre pos 2096 0.783 0.433 0.433 ns 2226 1.066 0.287 0.287 ns
Controle 9 ano pre pos 2096 -0.131 0.896 0.896 ns 2226 0.246 0.806 0.806 ns
Controle 6 ano 7 ano 1047 1.644 0.100 0.602 ns 1112 0.755 0.451 1.000 ns
Controle 6 ano 8 ano 1047 -2.433 0.015 0.091 ns 1112 -2.513 0.012 0.073 ns
Controle 6 ano 9 ano 1047 -5.074 0.000 0.000 **** 1112 -4.792 0.000 0.000 ****
Controle 7 ano 8 ano 1047 -3.939 0.000 0.001 *** 1112 -3.235 0.001 0.008 **
Controle 7 ano 9 ano 1047 -6.756 0.000 0.000 **** 1112 -5.654 0.000 0.000 ****
Controle 8 ano 9 ano 1047 -2.300 0.022 0.130 ns 1112 -1.937 0.053 0.318 ns
Experimental 6 ano pre pos 2096 1.773 0.076 0.076 ns 2226 2.498 0.013 0.013 *
Experimental 7 ano pre pos 2096 0.944 0.345 0.345 ns 2226 1.665 0.096 0.096 ns
Experimental 8 ano pre pos 2096 -0.962 0.336 0.336 ns 2226 -0.986 0.324 0.324 ns
Experimental 9 ano pre pos 2096 -1.132 0.258 0.258 ns 2226 -0.596 0.551 0.551 ns
Experimental 6 ano 7 ano 1047 -2.574 0.010 0.061 ns 1112 -2.506 0.012 0.074 ns
Experimental 6 ano 8 ano 1047 -6.514 0.000 0.000 **** 1112 -6.849 0.000 0.000 ****
Experimental 6 ano 9 ano 1047 -7.114 0.000 0.000 **** 1112 -6.535 0.000 0.000 ****
Experimental 7 ano 8 ano 1047 -4.242 0.000 0.000 *** 1112 -4.734 0.000 0.000 ****
Experimental 7 ano 9 ano 1047 -4.930 0.000 0.000 **** 1112 -4.449 0.000 0.000 ****
Experimental 8 ano 9 ano 1047 -0.688 0.492 1.000 ns 1112 0.281 0.779 1.000 ns
6 ano Controle Experimental 1047 -0.249 0.803 0.803 ns 1112 0.121 0.904 0.904 ns
7 ano Controle Experimental 1047 -4.500 0.000 0.000 **** 1112 -3.134 0.002 0.002 **
8 ano Controle Experimental 1047 -3.348 0.001 0.001 *** 1112 -3.222 0.001 0.001 **
9 ano Controle Experimental 1047 -1.694 0.091 0.091 ns 1112 -1.070 0.285 0.285 ns
Controle 1st quintile pre pos 2014 1.024 0.306 0.306 ns 2222 0.622 0.534 0.534 ns
Controle 2nd quintile pre pos 2014 3.105 0.002 0.002 ** 2222 2.906 0.004 0.004 **
Controle 3rd quintile pre pos 2014 4.794 0.000 0.000 **** 2222 5.243 0.000 0.000 ****
Controle 4th quintile pre pos 2014 3.089 0.002 0.002 ** 2222 3.992 0.000 0.000 ****
Controle 5th quintile pre pos 2014 2.704 0.007 0.007 ** 2222 3.512 0.000 0.000 ***
Controle 1st quintile 2nd quintile 1006 1.565 0.118 1.000 ns 1110 2.058 0.040 0.398 ns
Controle 1st quintile 3rd quintile 1006 2.102 0.036 0.358 ns 1110 3.034 0.002 0.025 *
Controle 1st quintile 4th quintile 1006 -0.118 0.906 1.000 ns 1110 1.355 0.176 1.000 ns
Controle 1st quintile 5th quintile 1006 -0.370 0.711 1.000 ns 1110 1.061 0.289 1.000 ns
Controle 2nd quintile 3rd quintile 1006 1.201 0.230 1.000 ns 1110 1.940 0.053 0.527 ns
Controle 2nd quintile 4th quintile 1006 -1.486 0.138 1.000 ns 1110 0.038 0.970 1.000 ns
Controle 2nd quintile 5th quintile 1006 -1.515 0.130 1.000 ns 1110 0.000 1.000 1.000 ns
Controle 3rd quintile 4th quintile 1006 -3.359 0.001 0.008 ** 1110 -2.364 0.018 0.183 ns
Controle 3rd quintile 5th quintile 1006 -3.197 0.001 0.014 * 1110 -1.923 0.055 0.548 ns
Controle 4th quintile 5th quintile 1006 -0.703 0.482 1.000 ns 1110 -0.057 0.954 1.000 ns
Experimental 1st quintile pre pos 2014 -1.416 0.157 0.157 ns 2222 -1.406 0.160 0.160 ns
Experimental 2nd quintile pre pos 2014 0.824 0.410 0.410 ns 2222 1.055 0.292 0.292 ns
Experimental 3rd quintile pre pos 2014 0.152 0.879 0.879 ns 2222 2.465 0.014 0.014 *
Experimental 4th quintile pre pos 2014 -0.569 0.569 0.569 ns 2222 1.583 0.114 0.114 ns
Experimental 5th quintile pre pos 2014 1.120 0.263 0.263 ns 2222 2.436 0.015 0.015 *
Experimental 1st quintile 2nd quintile 1006 0.882 0.378 1.000 ns 1110 1.705 0.088 0.884 ns
Experimental 1st quintile 3rd quintile 1006 -0.005 0.996 1.000 ns 1110 1.916 0.056 0.556 ns
Experimental 1st quintile 4th quintile 1006 -0.482 0.630 1.000 ns 1110 1.324 0.186 1.000 ns
Experimental 1st quintile 5th quintile 1006 -0.099 0.921 1.000 ns 1110 1.342 0.180 1.000 ns
Experimental 2nd quintile 3rd quintile 1006 -1.014 0.311 1.000 ns 1110 1.008 0.314 1.000 ns
Experimental 2nd quintile 4th quintile 1006 -1.503 0.133 1.000 ns 1110 0.462 0.644 1.000 ns
Experimental 2nd quintile 5th quintile 1006 -0.722 0.470 1.000 ns 1110 0.750 0.453 1.000 ns
Experimental 3rd quintile 4th quintile 1006 -0.991 0.322 1.000 ns 1110 -0.437 0.662 1.000 ns
Experimental 3rd quintile 5th quintile 1006 -0.183 0.855 1.000 ns 1110 0.248 0.804 1.000 ns
Experimental 4th quintile 5th quintile 1006 0.707 0.480 1.000 ns 1110 0.730 0.466 1.000 ns
1st quintile Controle Experimental 1006 -2.176 0.030 0.030 * 1110 -1.642 0.101 0.101 ns
2nd quintile Controle Experimental 1006 -2.614 0.009 0.009 ** 1110 -2.002 0.046 0.046 *
3rd quintile Controle Experimental 1006 -4.957 0.000 0.000 **** 1110 -3.647 0.000 0.000 ***
4th quintile Controle Experimental 1006 -3.161 0.002 0.002 ** 1110 -1.944 0.052 0.052 ns
5th quintile Controle Experimental 1006 -1.453 0.147 0.147 ns 1110 -0.974 0.330 0.330 ns

EMMS Table Comparison

df <- do.call(plyr::rbind.fill, lemms)
df[["N-N'"]] <- df[["N"]] - df[["N'"]]
df <- df[,c(names(lfatores)[names(lfatores) %in% colnames(df)],
            names(df)[!names(df) %in% names(lfatores)])]
grupo Sexo Zona Cor.Raca Serie score.tde.quintile N M (pre) SE (pre) M (unadj) SE (unadj) M (adj) SE (adj) conf.low conf.high N’ M (pre)’ SE (pre)’ M (unadj)’ SE (unadj)’ M (adj)’ SE (adj)’ conf.low’ conf.high’ N-N’
Controle 424 37.351 0.975 35.432 1.018 35.112 0.358 34.409 35.815 485 37.631 0.871 33.816 0.953 33.348 0.480 32.405 34.290 -61
Experimental 557 36.770 0.741 37.230 0.777 37.473 0.313 36.860 38.087 636 36.748 0.670 35.376 0.749 35.733 0.419 34.910 36.556 -79
Controle F 208 40.635 1.278 38.615 1.351 35.106 0.503 34.119 36.093 247 40.652 1.098 36.567 1.265 33.273 0.676 31.945 34.600 -39
Controle M 211 33.934 1.457 32.393 1.517 35.389 0.499 34.410 36.368 238 34.496 1.332 30.962 1.410 33.427 0.688 32.077 34.776 -27
Experimental F 276 39.623 1.053 40.272 1.123 37.744 0.436 36.889 38.600 319 39.408 0.941 38.138 1.086 36.008 0.594 34.842 37.173 -43
Experimental M 279 34.082 1.030 34.484 1.053 37.335 0.434 36.483 38.187 317 34.073 0.931 32.596 1.011 35.457 0.597 34.285 36.628 -38
Controle Rural 228 35.956 1.330 34.978 1.426 35.358 0.547 34.284 36.432 243 36.181 1.259 34.243 1.380 34.462 0.669 33.149 35.775 -15
Controle Urbana 101 39.743 1.823 31.099 2.044 27.763 0.824 26.145 29.380 109 39.468 1.727 29.835 1.951 26.939 1.001 24.975 28.903 -8
Experimental Rural 260 34.215 1.113 34.777 1.186 36.865 0.514 35.857 37.874 284 34.592 1.036 33.627 1.141 35.352 0.620 34.136 36.569 -24
Experimental Urbana 151 38.318 1.389 36.099 1.501 34.161 0.673 32.840 35.483 167 37.850 1.287 34.108 1.462 32.745 0.807 31.160 34.330 -16
Controle Branca 45 41.378 2.626 36.600 3.007 30.951 1.399 28.202 33.701 50 41.460 2.396 34.300 2.902 29.123 1.619 25.941 32.304 -5
Controle Indígena 11 42.000 5.215 45.636 6.009 39.393 2.820 33.851 44.935 11 42.000 5.215 45.636 6.009 39.964 3.439 33.206 46.722 0
Controle Parda 150 36.320 1.578 33.887 1.695 33.075 0.763 31.576 34.574 162 36.741 1.485 32.463 1.634 31.613 0.895 29.854 33.373 -12
Experimental Branca 58 34.517 2.651 34.586 2.746 35.499 1.226 33.088 37.909 61 34.623 2.526 33.311 2.713 34.404 1.459 31.537 37.270 -3
Experimental Indígena 15 27.067 4.406 30.333 4.574 38.371 2.420 33.616 43.126 15 27.067 4.406 30.333 4.574 38.355 2.952 32.554 44.155 0
Experimental Parda 175 33.851 1.343 33.514 1.460 35.064 0.707 33.674 36.453 190 34.253 1.272 32.447 1.412 33.879 0.828 32.253 35.505 -15
Controle 6 ano 126 28.444 1.766 24.619 1.828 32.399 0.756 30.915 33.882 134 28.993 1.687 24.231 1.741 31.423 0.889 29.678 33.168 -8
Controle 7 ano 127 34.024 1.632 28.000 1.773 30.663 0.744 29.204 32.123 141 34.716 1.503 28.362 1.739 30.495 0.856 28.815 32.175 -14
Controle 8 ano 86 42.430 1.785 40.337 1.933 35.292 0.906 33.514 37.069 89 42.584 1.737 39.787 1.897 34.967 1.081 32.846 37.087 -3
Controle 9 ano 116 46.612 1.611 46.914 1.537 38.034 0.790 36.483 39.584 121 46.950 1.557 46.397 1.495 37.718 0.940 35.874 39.562 -5
Experimental 6 ano 150 28.607 1.324 25.020 1.367 32.651 0.694 31.288 34.013 159 29.000 1.265 24.094 1.329 31.279 0.819 29.673 32.885 -9
Experimental 7 ano 170 36.294 1.168 34.500 1.284 35.081 0.642 33.822 36.340 187 36.540 1.089 33.524 1.254 34.046 0.743 32.588 35.503 -17
Experimental 8 ano 141 39.248 1.525 41.255 1.600 39.128 0.705 37.744 40.512 143 39.000 1.516 41.042 1.604 39.390 0.850 37.722 41.057 -2
Experimental 9 ano 140 43.143 1.354 45.514 1.288 39.815 0.713 38.416 41.215 147 43.204 1.307 44.422 1.331 39.054 0.845 37.397 40.711 -7
Controle 1st quintile 112 9.152 0.522 8.205 0.750 34.585 1.783 31.086 38.083 113 9.230 0.523 8.522 0.808 37.046 2.325 32.484 41.608 -1
Controle 2nd quintile 51 24.392 0.530 20.137 1.288 32.089 1.322 29.495 34.682 59 24.542 0.492 19.966 1.686 32.835 1.666 29.566 36.104 -8
Controle 3rd quintile 42 38.905 0.489 31.667 1.800 29.879 1.216 27.493 32.265 54 38.815 0.409 30.185 2.082 28.463 1.437 25.644 31.282 -12
Controle 4th quintile 125 47.544 0.223 44.840 0.850 34.874 0.932 33.046 36.702 135 47.370 0.222 43.215 0.978 32.746 1.190 30.411 35.080 -10
Controle 5th quintile 116 58.819 0.505 56.362 0.792 35.722 1.463 32.851 38.593 124 58.621 0.482 54.806 0.979 32.835 1.876 29.155 36.516 -8
Experimental 1st quintile 108 9.602 0.522 10.935 0.788 36.889 1.765 33.425 40.352 112 9.777 0.515 11.384 0.898 39.349 2.290 34.856 43.842 -4
Experimental 2nd quintile 101 25.842 0.347 25.040 0.976 35.619 1.016 33.625 37.613 117 25.855 0.316 24.675 1.128 36.203 1.292 33.668 38.737 -16
Experimental 3rd quintile 116 37.371 0.301 37.233 0.859 36.898 0.729 35.467 38.328 141 37.383 0.268 34.872 1.096 34.614 0.886 32.876 36.352 -25
Experimental 4th quintile 120 46.800 0.230 47.308 0.725 38.047 0.915 36.251 39.842 132 46.697 0.217 45.030 0.924 35.250 1.166 32.963 37.537 -12
Experimental 5th quintile 126 58.563 0.448 57.587 0.700 37.189 1.435 34.372 40.006 134 58.336 0.434 55.791 0.917 34.111 1.839 30.503 37.720 -8

CONSORT tables

Enrolment participants table

df <- do.call(rbind, lapply(c("Sexo","Zona","Cor.Raca","Serie","Idade"), FUN = function(x) {
  tdat2 <- tdat
  tdat2[[x]][is.na(tdat2[[x]])] <- "Not declared" 
  count_(group_by(tdat2, grupo), x)
}))
## Warning: `count_()` was deprecated in dplyr 0.7.0.
## ℹ Please use `count()` instead.
## ℹ See vignette('programming') for more help
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

## Warning: `count_()` was deprecated in dplyr 0.7.0.
## ℹ Please use `count()` instead.
## ℹ See vignette('programming') for more help
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

## Warning: `count_()` was deprecated in dplyr 0.7.0.
## ℹ Please use `count()` instead.
## ℹ See vignette('programming') for more help
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

## Warning: `count_()` was deprecated in dplyr 0.7.0.
## ℹ Please use `count()` instead.
## ℹ See vignette('programming') for more help
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

## Warning: `count_()` was deprecated in dplyr 0.7.0.
## ℹ Please use `count()` instead.
## ℹ See vignette('programming') for more help
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
grupo Sexo n Zona Cor.Raca Serie Idade
Controle F 376
Controle M 393
Experimental F 508
Experimental M 528
Controle 211 Not declared
Controle 351 Rural
Controle 207 Urbana
Experimental 321 Not declared
Experimental 426 Rural
Experimental 289 Urbana
Controle 77 Branca
Controle 21 Indígena
Controle 416 Not declared
Controle 254 Parda
Controle 1 Preta
Experimental 3 Amarela
Experimental 92 Branca
Experimental 22 Indígena
Experimental 619 Not declared
Experimental 299 Parda
Experimental 1 Preta
Controle 215 6 ano
Controle 214 7 ano
Controle 163 8 ano
Controle 177 9 ano
Experimental 250 6 ano
Experimental 283 7 ano
Experimental 258 8 ano
Experimental 245 9 ano
Controle 112 11
Controle 162 12
Controle 183 13
Controle 172 14
Controle 90 15
Controle 30 16
Controle 9 17
Controle 5 18
Controle 1 19
Controle 2 20
Controle 1 21
Controle 1 22
Controle 1 34
Experimental 1 1
Experimental 160 11
Experimental 215 12
Experimental 226 13
Experimental 220 14
Experimental 126 15
Experimental 52 16
Experimental 21 17
Experimental 7 18
Experimental 3 19
Experimental 1 20
Experimental 1 21
Experimental 1 23
Experimental 1 4
Experimental 1 9

Allocated participants table (after removed students with disabilities/special needs)

df <- do.call(rbind, lapply(c("Sexo","Zona","Cor.Raca","Serie","Idade"), FUN = function(x) {
  tdat2 <- gdat
  tdat2[[x]][is.na(tdat2[[x]])] <- "Not declared" 
  count_(group_by(tdat2, grupo), x)
}))
## Warning: `count_()` was deprecated in dplyr 0.7.0.
## ℹ Please use `count()` instead.
## ℹ See vignette('programming') for more help
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

## Warning: `count_()` was deprecated in dplyr 0.7.0.
## ℹ Please use `count()` instead.
## ℹ See vignette('programming') for more help
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

## Warning: `count_()` was deprecated in dplyr 0.7.0.
## ℹ Please use `count()` instead.
## ℹ See vignette('programming') for more help
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

## Warning: `count_()` was deprecated in dplyr 0.7.0.
## ℹ Please use `count()` instead.
## ℹ See vignette('programming') for more help
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

## Warning: `count_()` was deprecated in dplyr 0.7.0.
## ℹ Please use `count()` instead.
## ℹ See vignette('programming') for more help
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
grupo Sexo n Zona Cor.Raca Serie Idade
Controle F 367
Controle M 386
Experimental F 502
Experimental M 515
Controle 208 Not declared
Controle 343 Rural
Controle 202 Urbana
Experimental 315 Not declared
Experimental 421 Rural
Experimental 281 Urbana
Controle 76 Branca
Controle 20 Indígena
Controle 404 Not declared
Controle 252 Parda
Controle 1 Preta
Experimental 2 Amarela
Experimental 89 Branca
Experimental 22 Indígena
Experimental 610 Not declared
Experimental 293 Parda
Experimental 1 Preta
Controle 213 6 ano
Controle 209 7 ano
Controle 158 8 ano
Controle 173 9 ano
Experimental 245 6 ano
Experimental 277 7 ano
Experimental 253 8 ano
Experimental 242 9 ano
Controle 112 11
Controle 161 12
Controle 181 13
Controle 166 14
Controle 87 15
Controle 30 16
Controle 7 17
Controle 5 18
Controle 1 20
Controle 1 21
Controle 1 22
Controle 1 34
Experimental 1 1
Experimental 157 11
Experimental 214 12
Experimental 220 13
Experimental 218 14
Experimental 124 15
Experimental 50 16
Experimental 19 17
Experimental 6 18
Experimental 3 19
Experimental 1 20
Experimental 1 21
Experimental 1 23
Experimental 1 4
Experimental 1 9

Follow-up participants table (after removed students with does-not have pre and post-test)

df <- do.call(rbind, lapply(c("Sexo","Zona","Cor.Raca","Serie","Idade"), FUN = function(x) {
  tdat2 <- gdat[gdat$id %in% dat$id,]  
  tdat2[[x]][is.na(tdat2[[x]])] <- "Not declared" 
  count_(group_by(tdat2, grupo), x)
}))
## Warning: `count_()` was deprecated in dplyr 0.7.0.
## ℹ Please use `count()` instead.
## ℹ See vignette('programming') for more help
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

## Warning: `count_()` was deprecated in dplyr 0.7.0.
## ℹ Please use `count()` instead.
## ℹ See vignette('programming') for more help
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

## Warning: `count_()` was deprecated in dplyr 0.7.0.
## ℹ Please use `count()` instead.
## ℹ See vignette('programming') for more help
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

## Warning: `count_()` was deprecated in dplyr 0.7.0.
## ℹ Please use `count()` instead.
## ℹ See vignette('programming') for more help
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

## Warning: `count_()` was deprecated in dplyr 0.7.0.
## ℹ Please use `count()` instead.
## ℹ See vignette('programming') for more help
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
grupo Sexo n Zona Cor.Raca Serie Idade
Controle F 247
Controle M 238
Experimental F 319
Experimental M 317
Controle 133 Not declared
Controle 243 Rural
Controle 109 Urbana
Experimental 185 Not declared
Experimental 284 Rural
Experimental 167 Urbana
Controle 50 Branca
Controle 11 Indígena
Controle 262 Not declared
Controle 162 Parda
Experimental 1 Amarela
Experimental 61 Branca
Experimental 15 Indígena
Experimental 369 Not declared
Experimental 190 Parda
Controle 134 6 ano
Controle 141 7 ano
Controle 89 8 ano
Controle 121 9 ano
Experimental 159 6 ano
Experimental 187 7 ano
Experimental 143 8 ano
Experimental 147 9 ano
Controle 85 11
Controle 109 12
Controle 115 13
Controle 108 14
Controle 45 15
Controle 14 16
Controle 5 17
Controle 1 18
Controle 1 21
Controle 1 22
Controle 1 34
Experimental 1 1
Experimental 112 11
Experimental 153 12
Experimental 139 13
Experimental 132 14
Experimental 68 15
Experimental 21 16
Experimental 6 17
Experimental 1 18
Experimental 1 19
Experimental 1 21
Experimental 1 9

Participants table for Analysis (after exclude non-normal data)

df <- do.call(rbind, lapply(names(ldat), FUN = function(tname) {
  dat2 <- ldat[[tname]]
  data.frame(
    "For Analysis of ANCOVA" = tname,
    do.call(rbind, lapply(c("Sexo","Zona","Cor.Raca","Serie","Idade"), FUN = function(x) {
      tdat2 <- gdat[gdat$id %in% dat2$id,]  
      tdat2[[x]][is.na(tdat2[[x]])] <- "Not declared" 
      count_(group_by(tdat2, grupo), x)
  })))
}))  
## Warning: `count_()` was deprecated in dplyr 0.7.0.
## ℹ Please use `count()` instead.
## ℹ See vignette('programming') for more help
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

## Warning: `count_()` was deprecated in dplyr 0.7.0.
## ℹ Please use `count()` instead.
## ℹ See vignette('programming') for more help
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

## Warning: `count_()` was deprecated in dplyr 0.7.0.
## ℹ Please use `count()` instead.
## ℹ See vignette('programming') for more help
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

## Warning: `count_()` was deprecated in dplyr 0.7.0.
## ℹ Please use `count()` instead.
## ℹ See vignette('programming') for more help
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

## Warning: `count_()` was deprecated in dplyr 0.7.0.
## ℹ Please use `count()` instead.
## ℹ See vignette('programming') for more help
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

## Warning: `count_()` was deprecated in dplyr 0.7.0.
## ℹ Please use `count()` instead.
## ℹ See vignette('programming') for more help
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

## Warning: `count_()` was deprecated in dplyr 0.7.0.
## ℹ Please use `count()` instead.
## ℹ See vignette('programming') for more help
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

## Warning: `count_()` was deprecated in dplyr 0.7.0.
## ℹ Please use `count()` instead.
## ℹ See vignette('programming') for more help
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

## Warning: `count_()` was deprecated in dplyr 0.7.0.
## ℹ Please use `count()` instead.
## ℹ See vignette('programming') for more help
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

## Warning: `count_()` was deprecated in dplyr 0.7.0.
## ℹ Please use `count()` instead.
## ℹ See vignette('programming') for more help
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

## Warning: `count_()` was deprecated in dplyr 0.7.0.
## ℹ Please use `count()` instead.
## ℹ See vignette('programming') for more help
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

## Warning: `count_()` was deprecated in dplyr 0.7.0.
## ℹ Please use `count()` instead.
## ℹ See vignette('programming') for more help
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

## Warning: `count_()` was deprecated in dplyr 0.7.0.
## ℹ Please use `count()` instead.
## ℹ See vignette('programming') for more help
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

## Warning: `count_()` was deprecated in dplyr 0.7.0.
## ℹ Please use `count()` instead.
## ℹ See vignette('programming') for more help
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

## Warning: `count_()` was deprecated in dplyr 0.7.0.
## ℹ Please use `count()` instead.
## ℹ See vignette('programming') for more help
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

## Warning: `count_()` was deprecated in dplyr 0.7.0.
## ℹ Please use `count()` instead.
## ℹ See vignette('programming') for more help
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

## Warning: `count_()` was deprecated in dplyr 0.7.0.
## ℹ Please use `count()` instead.
## ℹ See vignette('programming') for more help
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

## Warning: `count_()` was deprecated in dplyr 0.7.0.
## ℹ Please use `count()` instead.
## ℹ See vignette('programming') for more help
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

## Warning: `count_()` was deprecated in dplyr 0.7.0.
## ℹ Please use `count()` instead.
## ℹ See vignette('programming') for more help
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

## Warning: `count_()` was deprecated in dplyr 0.7.0.
## ℹ Please use `count()` instead.
## ℹ See vignette('programming') for more help
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

## Warning: `count_()` was deprecated in dplyr 0.7.0.
## ℹ Please use `count()` instead.
## ℹ See vignette('programming') for more help
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

## Warning: `count_()` was deprecated in dplyr 0.7.0.
## ℹ Please use `count()` instead.
## ℹ See vignette('programming') for more help
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

## Warning: `count_()` was deprecated in dplyr 0.7.0.
## ℹ Please use `count()` instead.
## ℹ See vignette('programming') for more help
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

## Warning: `count_()` was deprecated in dplyr 0.7.0.
## ℹ Please use `count()` instead.
## ℹ See vignette('programming') for more help
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

## Warning: `count_()` was deprecated in dplyr 0.7.0.
## ℹ Please use `count()` instead.
## ℹ See vignette('programming') for more help
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

## Warning: `count_()` was deprecated in dplyr 0.7.0.
## ℹ Please use `count()` instead.
## ℹ See vignette('programming') for more help
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

## Warning: `count_()` was deprecated in dplyr 0.7.0.
## ℹ Please use `count()` instead.
## ℹ See vignette('programming') for more help
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

## Warning: `count_()` was deprecated in dplyr 0.7.0.
## ℹ Please use `count()` instead.
## ℹ See vignette('programming') for more help
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

## Warning: `count_()` was deprecated in dplyr 0.7.0.
## ℹ Please use `count()` instead.
## ℹ See vignette('programming') for more help
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

## Warning: `count_()` was deprecated in dplyr 0.7.0.
## ℹ Please use `count()` instead.
## ℹ See vignette('programming') for more help
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
For.Analysis.of.ANCOVA grupo Sexo n Zona Cor.Raca Serie Idade
grupo Controle F 207
grupo Controle M 217
grupo Experimental F 279
grupo Experimental M 278
grupo Controle 122 Not declared
grupo Controle 220 Rural
grupo Controle 82 Urbana
grupo Experimental 168 Not declared
grupo Experimental 250 Rural
grupo Experimental 139 Urbana
grupo Controle 40 Branca
grupo Controle 11 Indígena
grupo Controle 239 Not declared
grupo Controle 134 Parda
grupo Experimental 1 Amarela
grupo Experimental 53 Branca
grupo Experimental 14 Indígena
grupo Experimental 332 Not declared
grupo Experimental 157 Parda
grupo Controle 119 6 ano
grupo Controle 106 7 ano
grupo Controle 82 8 ano
grupo Controle 117 9 ano
grupo Experimental 126 6 ano
grupo Experimental 158 7 ano
grupo Experimental 137 8 ano
grupo Experimental 136 9 ano
grupo Controle 78 11
grupo Controle 79 12
grupo Controle 100 13
grupo Controle 104 14
grupo Controle 42 15
grupo Controle 12 16
grupo Controle 5 17
grupo Controle 1 18
grupo Controle 1 21
grupo Controle 1 22
grupo Controle 1 34
grupo Experimental 1 1
grupo Experimental 83 11
grupo Experimental 138 12
grupo Experimental 126 13
grupo Experimental 119 14
grupo Experimental 63 15
grupo Experimental 19 16
grupo Experimental 5 17
grupo Experimental 1 18
grupo Experimental 1 19
grupo Experimental 1 21
grupo:Sexo Controle F 208
grupo:Sexo Controle M 211
grupo:Sexo Experimental F 276
grupo:Sexo Experimental M 279
grupo:Sexo Controle 119 Not declared
grupo:Sexo Controle 219 Rural
grupo:Sexo Controle 81 Urbana
grupo:Sexo Experimental 164 Not declared
grupo:Sexo Experimental 248 Rural
grupo:Sexo Experimental 143 Urbana
grupo:Sexo Controle 39 Branca
grupo:Sexo Controle 11 Indígena
grupo:Sexo Controle 234 Not declared
grupo:Sexo Controle 135 Parda
grupo:Sexo Experimental 1 Amarela
grupo:Sexo Experimental 52 Branca
grupo:Sexo Experimental 15 Indígena
grupo:Sexo Experimental 329 Not declared
grupo:Sexo Experimental 158 Parda
grupo:Sexo Controle 117 6 ano
grupo:Sexo Controle 106 7 ano
grupo:Sexo Controle 82 8 ano
grupo:Sexo Controle 114 9 ano
grupo:Sexo Experimental 123 6 ano
grupo:Sexo Experimental 158 7 ano
grupo:Sexo Experimental 137 8 ano
grupo:Sexo Experimental 137 9 ano
grupo:Sexo Controle 77 11
grupo:Sexo Controle 80 12
grupo:Sexo Controle 99 13
grupo:Sexo Controle 102 14
grupo:Sexo Controle 41 15
grupo:Sexo Controle 11 16
grupo:Sexo Controle 5 17
grupo:Sexo Controle 1 18
grupo:Sexo Controle 1 21
grupo:Sexo Controle 1 22
grupo:Sexo Controle 1 34
grupo:Sexo Experimental 1 1
grupo:Sexo Experimental 80 11
grupo:Sexo Experimental 138 12
grupo:Sexo Experimental 125 13
grupo:Sexo Experimental 120 14
grupo:Sexo Experimental 64 15
grupo:Sexo Experimental 19 16
grupo:Sexo Experimental 5 17
grupo:Sexo Experimental 1 18
grupo:Sexo Experimental 1 19
grupo:Sexo Experimental 1 21
grupo:Zona Controle F 164
grupo:Zona Controle M 165
grupo:Zona Experimental F 201
grupo:Zona Experimental M 210
grupo:Zona Controle 228 Rural
grupo:Zona Controle 101 Urbana
grupo:Zona Experimental 260 Rural
grupo:Zona Experimental 151 Urbana
grupo:Zona Controle 37 Branca
grupo:Zona Controle 11 Indígena
grupo:Zona Controle 167 Not declared
grupo:Zona Controle 114 Parda
grupo:Zona Experimental 46 Branca
grupo:Zona Experimental 15 Indígena
grupo:Zona Experimental 205 Not declared
grupo:Zona Experimental 145 Parda
grupo:Zona Controle 110 6 ano
grupo:Zona Controle 110 7 ano
grupo:Zona Controle 49 8 ano
grupo:Zona Controle 60 9 ano
grupo:Zona Experimental 110 6 ano
grupo:Zona Experimental 145 7 ano
grupo:Zona Experimental 73 8 ano
grupo:Zona Experimental 83 9 ano
grupo:Zona Controle 75 11
grupo:Zona Controle 85 12
grupo:Zona Controle 73 13
grupo:Zona Controle 56 14
grupo:Zona Controle 25 15
grupo:Zona Controle 9 16
grupo:Zona Controle 2 17
grupo:Zona Controle 1 18
grupo:Zona Controle 1 21
grupo:Zona Controle 1 22
grupo:Zona Controle 1 34
grupo:Zona Experimental 1 1
grupo:Zona Experimental 73 11
grupo:Zona Experimental 132 12
grupo:Zona Experimental 80 13
grupo:Zona Experimental 64 14
grupo:Zona Experimental 43 15
grupo:Zona Experimental 12 16
grupo:Zona Experimental 4 17
grupo:Zona Experimental 1 19
grupo:Zona Experimental 1 21
grupo:Cor.Raca Controle F 107
grupo:Cor.Raca Controle M 99
grupo:Cor.Raca Experimental F 121
grupo:Cor.Raca Experimental M 127
grupo:Cor.Raca Controle 48 Not declared
grupo:Cor.Raca Controle 127 Rural
grupo:Cor.Raca Controle 31 Urbana
grupo:Cor.Raca Experimental 34 Not declared
grupo:Cor.Raca Experimental 173 Rural
grupo:Cor.Raca Experimental 41 Urbana
grupo:Cor.Raca Controle 45 Branca
grupo:Cor.Raca Controle 11 Indígena
grupo:Cor.Raca Controle 150 Parda
grupo:Cor.Raca Experimental 58 Branca
grupo:Cor.Raca Experimental 15 Indígena
grupo:Cor.Raca Experimental 175 Parda
grupo:Cor.Raca Controle 56 6 ano
grupo:Cor.Raca Controle 54 7 ano
grupo:Cor.Raca Controle 37 8 ano
grupo:Cor.Raca Controle 59 9 ano
grupo:Cor.Raca Experimental 71 6 ano
grupo:Cor.Raca Experimental 53 7 ano
grupo:Cor.Raca Experimental 69 8 ano
grupo:Cor.Raca Experimental 55 9 ano
grupo:Cor.Raca Controle 41 11
grupo:Cor.Raca Controle 38 12
grupo:Cor.Raca Controle 52 13
grupo:Cor.Raca Controle 50 14
grupo:Cor.Raca Controle 14 15
grupo:Cor.Raca Controle 7 16
grupo:Cor.Raca Controle 3 17
grupo:Cor.Raca Controle 1 21
grupo:Cor.Raca Experimental 52 11
grupo:Cor.Raca Experimental 47 12
grupo:Cor.Raca Experimental 52 13
grupo:Cor.Raca Experimental 50 14
grupo:Cor.Raca Experimental 30 15
grupo:Cor.Raca Experimental 14 16
grupo:Cor.Raca Experimental 2 17
grupo:Cor.Raca Experimental 1 21
grupo:Serie Controle F 230
grupo:Serie Controle M 225
grupo:Serie Experimental F 301
grupo:Serie Experimental M 300
grupo:Serie Controle 129 Not declared
grupo:Serie Controle 227 Rural
grupo:Serie Controle 99 Urbana
grupo:Serie Experimental 176 Not declared
grupo:Serie Experimental 270 Rural
grupo:Serie Experimental 155 Urbana
grupo:Serie Controle 48 Branca
grupo:Serie Controle 11 Indígena
grupo:Serie Controle 249 Not declared
grupo:Serie Controle 147 Parda
grupo:Serie Experimental 1 Amarela
grupo:Serie Experimental 57 Branca
grupo:Serie Experimental 15 Indígena
grupo:Serie Experimental 353 Not declared
grupo:Serie Experimental 175 Parda
grupo:Serie Controle 126 6 ano
grupo:Serie Controle 127 7 ano
grupo:Serie Controle 86 8 ano
grupo:Serie Controle 116 9 ano
grupo:Serie Experimental 150 6 ano
grupo:Serie Experimental 170 7 ano
grupo:Serie Experimental 141 8 ano
grupo:Serie Experimental 140 9 ano
grupo:Serie Controle 81 11
grupo:Serie Controle 99 12
grupo:Serie Controle 106 13
grupo:Serie Controle 104 14
grupo:Serie Controle 42 15
grupo:Serie Controle 14 16
grupo:Serie Controle 5 17
grupo:Serie Controle 1 18
grupo:Serie Controle 1 21
grupo:Serie Controle 1 22
grupo:Serie Controle 1 34
grupo:Serie Experimental 1 1
grupo:Serie Experimental 104 11
grupo:Serie Experimental 147 12
grupo:Serie Experimental 133 13
grupo:Serie Experimental 122 14
grupo:Serie Experimental 65 15
grupo:Serie Experimental 21 16
grupo:Serie Experimental 5 17
grupo:Serie Experimental 1 18
grupo:Serie Experimental 1 19
grupo:Serie Experimental 1 21
grupo:score.tde.quintile Controle F 223
grupo:score.tde.quintile Controle M 223
grupo:score.tde.quintile Experimental F 285
grupo:score.tde.quintile Experimental M 286
grupo:score.tde.quintile Controle 126 Not declared
grupo:score.tde.quintile Controle 227 Rural
grupo:score.tde.quintile Controle 93 Urbana
grupo:score.tde.quintile Experimental 171 Not declared
grupo:score.tde.quintile Experimental 256 Rural
grupo:score.tde.quintile Experimental 144 Urbana
grupo:score.tde.quintile Controle 44 Branca
grupo:score.tde.quintile Controle 10 Indígena
grupo:score.tde.quintile Controle 247 Not declared
grupo:score.tde.quintile Controle 145 Parda
grupo:score.tde.quintile Experimental 1 Amarela
grupo:score.tde.quintile Experimental 56 Branca
grupo:score.tde.quintile Experimental 15 Indígena
grupo:score.tde.quintile Experimental 336 Not declared
grupo:score.tde.quintile Experimental 163 Parda
grupo:score.tde.quintile Controle 123 6 ano
grupo:score.tde.quintile Controle 118 7 ano
grupo:score.tde.quintile Controle 87 8 ano
grupo:score.tde.quintile Controle 118 9 ano
grupo:score.tde.quintile Experimental 138 6 ano
grupo:score.tde.quintile Experimental 160 7 ano
grupo:score.tde.quintile Experimental 136 8 ano
grupo:score.tde.quintile Experimental 137 9 ano
grupo:score.tde.quintile Controle 79 11
grupo:score.tde.quintile Controle 92 12
grupo:score.tde.quintile Controle 104 13
grupo:score.tde.quintile Controle 106 14
grupo:score.tde.quintile Controle 43 15
grupo:score.tde.quintile Controle 13 16
grupo:score.tde.quintile Controle 5 17
grupo:score.tde.quintile Controle 1 18
grupo:score.tde.quintile Controle 1 21
grupo:score.tde.quintile Controle 1 22
grupo:score.tde.quintile Controle 1 34
grupo:score.tde.quintile Experimental 1 1
grupo:score.tde.quintile Experimental 92 11
grupo:score.tde.quintile Experimental 142 12
grupo:score.tde.quintile Experimental 126 13
grupo:score.tde.quintile Experimental 118 14
grupo:score.tde.quintile Experimental 64 15
grupo:score.tde.quintile Experimental 20 16
grupo:score.tde.quintile Experimental 5 17
grupo:score.tde.quintile Experimental 1 18
grupo:score.tde.quintile Experimental 1 19
grupo:score.tde.quintile Experimental 1 21